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SUMMARY 

An  integrated  modeling  and  simulation  program  has  been  conducted  to  substantially 
improve  the  fundamental  knowledge  of  supercritical  combustion  of  liquid  propellants  under 
conditions  representative  of  contemporary  rocket  engines.  Both  shear  and  swirl  co-axial  injectors 
were  considered.  The  formulation  was  based  on  the  complete  conservation  equations  for  a  multi- 
component  chemically  reacting  mixture  in  three  dimensions.  In  addition,  general-fluid 
thermodynamics  and  transport  theories  were  incorporated  to  allow  for  a  unified  treatment  of  fluid 
properties  over  the  entire  range  of  thermodynamic  states.  Turbulence  closure  was  achieved  by 
means  of  the  large-eddy-simulation  (LES)  technique.  Start-of-the-art  closure  schemes  for  subgrid- 
scale  dynamics  and  turbulence/chemistry  interactions  were  implemented.  Special  attention  was 
given  to  the  fluid  behavior  in  the  two-phase  and  transcritical  regimes  in  which  rapid  property 
variations  occur.  Furthermore,  an  efficient  numerical  framework  utilizing  contemporary  computer 
software  and  hardware  technologies  were  employed,  such  that  sweeping  calculations  were 
performed  within  a  realistic  time  frame. 

The  work  contained  the  following  four  tasks  in  a  hierarchical  manner: 

4.  development  of  an  efficient  numerical  algorithm  for  treating  general-fluid  mixtures  with 
chemical  reactions; 

5.  cryogenic  fluid  injection  and  mixing  at  supercritical  conditions; 

6.  shear  co-axial  injection  of  cryogenic  fluids  with  acoustic  excitations;  and 

7.  flow  and  flame  dynamics  of  shear  co-axial  injectors  with  liquid  oxygen  (LOX)  and 
methane 

Various  underlying  physiochemical  mechanisms  associated  with  co-axial  injector  dynamics 
were  studied  in  detail.  These  included  flow  evolution,  flame  stabilization  and  spreading,  heat 
transfer,  and  acoustic  response.  The  effects  of  design  attributes  (e.g.,  injection  port  size  and 
location,  center  post  recess  distance,  etc.)  and  operating  conditions  (e.g.,  chamber  pressure, 
velocity,  and  temperature,  swirl  strength,  etc.)  on  injector  characteristics  were  assessed.  Results 
have  not  only  enhanced  the  basic  understanding  of  the  subject  problem,  but  also  provided  a 
quantitative  basis  for  identifying  and  prioritizing  the  key  design  parameters  and  flow  variables  that 
exert  dominant  influences  on  the  injector  behavior  in  different  environments. 

Results  from  this  research  project  have  led  to  the  following  publications. 

1 .  “Dynamics  of  Simplex  Swirl  Injectors  for  Cryogenic  Propellants  at  Supercritical 
Conditions,”  by  N.  Zong  and  V.  Yang,  AIAA  Paper  2004-1332,  presented  at  the 
AIAA  42nd  Aerospace  Sciences  Meeting  and  Exhibit,  January  2004,  also 
submitted  to  Physics  of  Fluids,  2007. 

2.  “A  Numerical  Study  of  Cryogenic  Fluid  Injection  and  Mixing  under 
Supercritical  Conditions,”  by  N.  Zong,  H.  Meng,  S.Y.  Hsieh,  and  V.  Yang, 

Physics  of  Fluids,  Vol.  16,  2004,  pp.  4248-4261. 
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Task  1 

Development  of  an  Efficient  Numerical  Algorithm  for  Treating 
General-Fluid  Mixtures  with  Chemical  Reactions 

The  unified  treatment  of  general  fluid  thermodynamics  developed  by  Meng  and  Yang  for 
supercritical  fluid  transport  and  combustion  has  been  improved  and  incorporated  into  a  highly 
efficient  preconditioning  scheme  using  primitive  pressure-temperature  type  variables.  The  major 
advantageous  of  this  newly  developed  algorithm  include  1).  all  the  numerical  relationships, 
including  the  preconditioning  matrix,  Jacobian  matrices,  and  eigenvalues,  are  derived  based 
fundamental  thermodynamics  theories  with  the  concepts  of  partial-mass  and  partial-density 
properties.  This  renders  a  self-consistent  and  robust  algorithm  valid  for  fluid  flows  at  all  speeds 
and  at  all  thermodynamic  states;  2).  since  pressure,  temperature,  and  mass  fraction  are  selected  as 
the  primary  dependent  variables  of  the  governing  system  and  solved  directly,  the  cumbersome 
iterative  solution  of  real-fluid  equations  of  state  is  eliminated;  and  3).  the  original  governing  system 
can  be  fully  recovered  for  high-speed  flows  because  of  the  use  of  a  universal  and  general 
preconditioning  matrix. 

The  work  has  remedied  several  deficiencies  in  the  previous  studies.  Most  importantly, 
computational  efficiency  has  been  significantly  improved  because  the  tedious  iterative  solution  of  a 
real-fluid  equation  of  state  is  avoided  through  a  careful  selection  of  the  primary  dependent  variables 
of  the  governing  system.  In  addition,  load  balance  can  be  easily  achieved  on  a  distributed 
computing  facility.  The  resultant  scheme  is  highly  efficient  and  suitable  for  parallel 
implementation.  This  feature  is  extremely  attractive  because  of  the  large-scale  numerical 
calculations  required  for  three-dimensional  simulations  of  injector  flow  and  flame  dynamics.  A 
comprehensive  description  of  this  new  approach  is  given  in  the  following  paper  which  has  been 
submitted  for  publication  in  the  International  Journal  of  Computational  Fluid  Dynamics. 

“An  Efficient  Preconditioning  Scheme  for  General  Fluid  Mixtures  using  Primitive 
Pressure-Temperature  Variables,”  by  Zong,  N.,  and  Yang,  V.,  submitted  for 
publication  in  International  Journal  of  Computational  Fluid  Dynamics,  February 
2007. 

1.1  Introduction 

The  development  of  an  efficient  numerical  algorithm  capable  of  handling  fluid  flows  over  a 
broad  range  of  thermodynamic  states  is  of  particular  importance  for  many  scientific  and  engineering 
problems  involving  fluid  phase  transition  from  a  subcritical  to  a  supercritical  state.  Notable 
examples  include  fluid  transport,  heat  transfer,  material  processing,  and  high-pressure  combustion 
for  propulsion  and  power-generation  applications.  Thermodynamic  non-idealities  and  transport 
anomalies  often  occur  under  these  conditions,  especially  during  the  transition  through  the 
transcritical  regime  [1].  Thus,  treating  fluid-state  transition  and  thermophysical-property  variations 
in  a  manner  consistent  with  the  intrinsic  characteristics  of  a  numerical  algorithm  is  critical  to 
achieving  numerical  efficiency  and  robustness. 

Several  attempts  have  been  made  to  treat  fluid  flows  at  different  thermodynamic  states. 
Merkel  et  al.  [2]  extended  a  preconditioning  scheme  for  ideal  gases  [3-8]  to  accommodate  arbitrary 
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equations  of  state.  The  resultant  scheme,  along  with  the  use  of  the  Soave-Redlich-Kwong  equation 
of  state,  has  been  employed  to  simulate  supercritical  hydrogen  flow  through  a  two-dimensional 
cascade  and  a  heated  rectangular  duct.  Although  encouraging  results  were  obtained  for  the  specific 
problems  considered,  no  information  was  given  about  the  evaluation  of  thermophysical  properties. 
The  relationship  between  numerical  properties  and  general  fluid  thermodynamics  was  not 
addressed. 

Edwards  et  al.  [9]  considered  real-fluid  flows  with  liquid- vapor  phase  transition  using  a  low- 
diffusion  flux-splitting  scheme  [10].  The  fluid  p-v-T  behavior  in  the  liquid  phase  was  modeled  with 
the  Sanchez-Lacombe  equation  of  state.  A  homogeneous  vapor-liquid  phase  equilibrium  model 
was  employed  to  study  the  vapor  cavitation  during  liquid  carbon  dioxide  expansion  through  a 
sharp-orifice  nozzle.  Numerical  experiments  demonstrated  the  effectiveness  of  the  method  in 
capturing  several  important  two-phase  flow  phenomena,  such  as  cavitation  bubbles  and  vapor- 
liquid  condensation  shocks.  The  scheme  was  later  extended  to  investigate  gas-solid  two-phase 
flows  in  fluidized  beds  under  different  flow  conditions  [11]. 

Oefelein  and  Yang  [12]  extended  the  preconditioning  methodology  described  in  [13]  to 
handle  multi-component,  dense-fluid  mixtures  with  finite-rate  chemical  kinetics  at  supercritical 
pressures.  Temporally  accurate  solutions  were  obtained  through  the  implementation  of  a  dual-time- 
stepping  integration  technique,  in  which  the  primitive  variables  of  pressure,  velocity  components, 
temperature,  and  mass  fraction,  were  chosen  as  the  dependent  variables  in  pseudo-time.  The  Soave- 
Redlich-Kwong  equation  of  state  was  employed  to  predict  the  fluid  volumetric  behavior  except  in 
regions  near  the  critical  point,  for  which  the  Benedict-Webb-Rubin  equation  of  state  was  used. 
Thermodynamics  properties,  such  as  the  enthalpy,  Gibbs  energy,  and  specific  heat,  were  obtained 
as  explicit  functions  of  temperature  and  pressure  by  using  Maxwell’s  relations  to  derive 
thermodynamic  departure  functions  [1].  Transport  properties  were  estimated  using  an  extended 
corresponding-state  principle.  The  combined  theoretical/numerical  framework  was  applied  to  study 
the  mixing  and  combustion  of  cryogenic  propellants  in  supercritical  environments  [12,14].  The 
properties  of  fluid  mixtures  (i.e.,  enthalpy,  and  internal  energy)  in  their  analysis,  however,  were 
evaluated  with  the  formulas  for  an  ideal-gas  mixture,  taking  as  the  mass-weighted  sum  of  the  values 
of  constituent  species  without  considering  their  interactions.  Since  the  coupling  among  molecules 
of  different  species  may  be  significant  in  real  fluids,  such  a  simplification  leads  to  inconsistency 
with  real-fluid  equations  of  state  and  other  property-evaluation  schemes. 

To  remedy  this  deficiency,  Meng  and  Yang  [15]  developed  a  unified  numerical  treatment  of 
general  fluid  thermodynamics.  The  concepts  of  partial  mass  and  partial  density  were  introduced  to 
evaluate  the  properties  of  real-fluid  mixtures.  All  of  the  thermodynamic  properties  and  numerical 
relations,  including  the  Jacobian  matrices  and  eigenvalues,  were  derived  directly  from  fundamental 
thermodynamics  theories,  rendering  a  self-consistent  and  robust  algorithm  valid  over  the  entire 
range  of  fluid  states.  Furthermore,  full  account  is  taken  of  transport  property  variations  as  functions 
of  fluid  density,  temperature,  and  composition.  The  resultant  numerical  properties  and  property 
evaluation  routines  were  incorporated  into  a  preconditioning  scheme  [13]  to  handle  fluid  flows  at 
all  speeds.  The  overall  approach  is  general  and  can  accommodate  any  type  of  equation  of  state. 

In  spite  of  its  effectiveness  in  treating  real-fluid  flows  with  thermodynamic  non-idealities 
and  transport  anomalies  in  the  transcritical  regime  [15,  16],  the  aforementioned  formulation  can  not, 
in  general,  be  solved  in  a  non-iterative  manner  because  the  specific  enthalpy  is  used  as  one  of  the 
dependent  variables  in  the  formulation  of  the  pseudo-time  derivatives.  Extensive  iterations  at  each 
time  step  and  grid  point  are  required  to  determine  the  fluid  temperature  from  the  specific  enthalpy. 
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Although  such  iterations  are  straightforward  to  accomplish  for  an  ideal-gas  mixture,  the  procedures 
become  quite  complicated  and  time-consuming  for  real-fluid  mixtures  due  to  the  complicated  form 
of  the  equation  of  state  and  property  evaluation  schemes.  The  iterative  solution  of  real-fluid 
equation  of  state  not  only  significantly  increases  the  computational  cost,  but  also  impairs  the 
scalability  of  a  parallelized  code. 

The  purpose  of  the  present  work  is  to  circumvent  this  limitation  by  replacing  specific 
enthalpy  with  temperature  as  a  primary  dependent  variable  in  the  numerical  formulation.  Such  an 
employment  of  the  pressure-temperature  type  of  primitive  variables  has  been  made  in  many 
existing  approaches  [3-8]  and  offers  the  following  advantageous.  First,  laborious  iterations  in 
calculating  temperature  from  specific  enthalpy  or  internal  energy  are  avoid,  rendering  a  highly 
efficient  algorithm  for  real-fluid  mixtures.  Second,  load  balance  is  easier  to  achieve  on  a  distributed 
computing  facility  since  no  iterative  solution  of  the  equation  of  state  is  required.  The  computation 
burden  at  each  spatial  grid  point  is  the  same.  Third,  computation  of  numerical  Jacobian  matrices  is 
simplified,  especially  within  the  context  of  the  general  fluid  thermodynamics.  In  addition,  the 
preconditioning  matrix  in  this  study  is  derived  by  analyzing  the  eigenvalues  of  the  governing 
system.  All  the  off-diagonal  terms  in  the  Jacobian  matrix  relating  the  conservative  to 
preconditioning  variables  are  retained.  The  original  governing  equations  can  thus  be  fully 
recovered  for  high-speed  flows  in  a  steady-state  simulation.  The  scheme  is  further  optimized  to 
facilitate  parallel  computation  by  treating  the  pseudo-time  and  spatial  derivatives  using  a  fourth- 
order  Runge-Kutta  (RK4)  scheme  [17]  and  an  explicit  fourth-order  flux-differencing  scheme  [18], 
respectively.  Because  the  time  advancement  is  fully  explicit  in  the  pseudo-time  space,  the 
approach  not  only  guarantees  high-order  numerical  accuracy  in  both  time  and  space,  but  also 
greatly  simplifies  the  implementation  of  a  parallelized  code. 

The  remainder  of  this  part  is  organized  as  follows.  In  Section  1.2,  we  briefly  summarize  the 
theoretical  framework,  including  the  governing  equations  and  basic  thermodynamics  theories  for 
general  fluid  mixtures.  Fundamental  thermodynamic  and  numerical  attributes  are  established. 
Section  1.3  deals  with  the  numerical  implementation  of  the  present  scheme.  A  stability  analysis  is 
carried  out  in  Section  1.4  to  characterize  the  stability  and  convergence  behavior  of  the  algorithm. 
Section  1.5  presents  results  from  selected  numerical  experiments,  along  with  an  assessment  of  the 
computational  efficiency.  Finally,  a  brief  summary  concludes  the  work. 


1.2  Theoretical  Formulation 

To  facilitate  discussion,  we  consider  the  two-dimensional  conservation  equations  of  mass, 
momentum,  energy,  and  species  transport  for  a  chemically  reacting  system  of  N  species. 
Following  the  approach  detailed  in  [4,6]  for  developing  a  preconditioning  scheme  for  fluid  flows  at 

all  speeds,  a  pseudo-time  derivative  of  the  form  TdQ/ dr  is  added  to  the  conservation  equations, 


r3j2  |  ae  |  9(£-£,)  |  3(f-F„)  . 

dr  dt  dx  dx 

where  T  represents  the  preconditioning  matrix  and  r  the  pseudo-time, 
vector,  Q ,  is  defined  as, 


(1) 

The  conservative  variable 


Q  =  yslp  pu  pv  pe -  pY,f . 


(2) 
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The  exponent  8  =  0  or  1  corresponds  to  a  planar  two-dimensional  or  an  axisymmetric  case, 
respectively.  Standard  notation  in  fluid  mechanics  is  used  here,  with  p,  ( u,v ),  et,  and  Y. 

denoting  the  density,  velocity  components,  specific  total  energy,  and  mass  fraction  of  species  i, 
respectively.  Explicit  expressions  of  the  convective  flux  vectors,  E  and  F ,  and  the  diffusion  flux 
vectors,  Ev ,  and  Fv  ,  are  given  in  [13-15].  The  source  term  in  Eq.  (1),  S ,  arises  from  chemical 
reactions  or  axisymmetric  geometry  [13]. 

Because  the  pseudo-time  derivative  in  Eq.  (1)  vanishes  at  convergence,  a  certain  amount  of 
liberty  can  be  taken  in  the  selection  of  the  pseudo-time  variables  [8].  Different  options  have  been 
proposed  for  the  primary  dependent  variables  in  the  preconditioning  scheme  for  ideal  gases.  Those 
include  conservative  and  several  sets  of  primitive  (e.g.,  pressure-entropy,  pressure-enthalpy,  and 
pressure-temperature)  variables  [8].  The  influences  of  selected  preconditioning  variables  on  the 
numerical  convergence  and  accuracy  for  low-Mach  number  aerodynamic  problems  were  examined 
by  Turkel  and  colleagues  [19,20].  Results  indicated  that  the  pressure-temperature  type  of  primitive 
variables,  which  had  long  been  a  favorable  choice  by  many  researchers,  led  to  fast  convergence  for 
simulating  steady-state  flows.  The  present  study  follows  the  same  approach,  and  the  pseudo-time 

variable  vector,  Q ,  is  thus  defined  as 

Q  =  [pg,u,v,T,Yi].  (3) 

The  use  of  the  gauge  pressure,  pg ,  taken  as  the  static  pressure  subtracted  by  a  reference 
pressure  [4,13],  is  crucial  for  two  reasons.  First,  it  allows  the  vectors  E,  F ,  F  ,  Fv,  and  S  to  be 
expressed  as  unique  functions  of  Q  [13].  Second,  pg  is  much  more  sensitive  to  the  solution  than 

the  density  at  low  Mach  numbers,  and  thus  provides  a  stronger  velocity-pressure  coupling  in  the 
momentum  balance  [6].  The  selection  of  the  velocity  components  (w,  v),  temperature,  T,  and 
species  mass  fraction,  Yn  helps  maintain  the  direction  of  the  diffusion  process,  simplifies  the 
structures  of  viscous  vectors,  and  reduces  the  computational  complexity  [5].  In  contrast  to  the 
schemes  with  the  pressure-enthalpy  type  of  primitive  variables,  all  the  thermodynamic  properties  of 
pressure,  temperature,  and  mass  fraction  are  directly  solved,  thereby  eliminating  tedious  iterative 
solutions  of  the  equations  of  state. 

1.2.1  Thermodynamic  Relations 

To  develop  a  robust  numerical  scheme  in  a  manner  consistent  with  real-fluid 
thermodynamics,  all  the  numerical  properties  and  relations,  including  the  preconditioning  matrix, 
Jacobian  matrices,  and  system  eigenvalues,  must  be  derived  in  accordance  with  fundamental 
thermodynamics  theories  described  in  Meng  and  Yang  [15].  The  procedure  should  conform  to  the 
concepts  of  partial  mass  and  partial  density  properties  for  fluid  mixtures,  along  with  the  use  of  an 
appropriate  equation  of  state. 

Four  thermodynamic  relationships  beyond  those  established  in  [15]  are  needed  to  derive  the 
numerical  properties.  Those  equations  define  density,  enthalpy,  and  internal  energy  as  functions  of 
such  dependent  variables  as  pressure,  temperature,  and  mass  fraction,  and  consequently  can  be  used 

to  determine  the  Jacobian  matrices  T  =  dQ/dQ,  A  =  dE/dQ  ,  B  =  dF/dQ,  and  D  =  dS/dQ  (see 
Appendix  A),  as  well  as  the  associated  eigen-properties. 


6 


Pressure  as  function  of  density,  temperature,  and  mass  fraction 

The  first  relationship  expresses  pressure  as  a  function  of  temperature,  density,  and  mass 
fractions.  For  a  general  fluid  mixture  consisting  of  N  species,  each  intensive  thermodynamic 
property  can  be  determined  by  the  other  N  + 1  properties  of  the  mixture.  This  leads  to  the 
following  expression: 

p  =  p(T,pi)  (4) 


where  i  =  \,---,N.  The  differential  form  of  Eq.  (4)  is 


(5) 


We  may  change  the  summation  on  the  right  hand  side  of  Eq.  (5)  from  1  through  N  to  1 
through  N to  obtain  the  following  relation: 


dp  =  A  dT  +  V  [(— )r  -  (-^-)r 
y  dT  p'  tf  dp/T  p'«  dpN  'P"N 


1=1 


VP, 


dp 


(6) 


Since  pt  =  pYn 

dpt  =  Yidp  +  pdYi 

substitution  of  Eq.  (7)  into  Eq.  (6)  yields 
dp  =  Aj-dT  +  AydYi  +  Apd  p 


where 


A  =  A 

Vri; 


(7) 

(8) 


(9a) 

(9b) 

(9c) 


It  has  been  shown  in  Ref.  [15]  that  the  speed  of  sound  for  a  general  fluid  mixture  can  be 
expressed  as 


-yA 

cv(dp)T*  7  p 


(10) 


7 


where  Cp  and  Cv  are  constant-pressure  and  constant-volume  specific  heats,  respectively,  and  y  is 
the  ratio  of  specific  heat. 

Internal  energy  as  function  of  temperature,  pressure,  and  mass  fraction 

The  second  relationship  expresses  internal  energy  as  a  function  of  pressure,  density,  and 
mass  fractions.  We  begin  with  the  following  functional  relationship: 


pe  =  pe(T, pi) 


(11) 


where  i  =  l, ■■■, N  ,  and  e  the  specific  internal  energy.  The  differential  form  of  Eq.  (11)  can  be 
written  as 


dpe  =  p 


f  de  s 

Jrf, 


( 


dpe 


N 

dT  +  Y 

pi 


dpi 


(12) 


'T  ,p  j^i 


The  derivative  in  the  first  term  on  the  right  hand  side  is  the  constant- volume  heat  capacity, 
Cv ,  and  that  in  the  second  term  is  the  partial-density  internal  energy,  e,  [15], 


e. 


dpe 


V BP 


Equation  (12)  thus  becomes 

N 

dpe  =  pCvdT  +  Y^^dPi 
1=1 

Substitution  of  Eq.  (7)  into  Eq.  (14)  gives 

dpe  =  pCvdT  +  YjeipdYi  +J^elYidp 


i=i 


Since  dpe  =  pde  +  edp ,  a  simple  manipulation  of  Eq.  (15)  results  in 

de  =  CvdT  +  §(«,  - eN )dYi  +-(fj Ye,  - e)dp 
i=l  P  i=l 

Substitution  of  Eq.  (8)  into  the  above  equation  leads  to 

N- 1 

de  =  BrdT  +  Bpdp  +  ^  Bv  dY, 

i=i 


(13) 


(14) 


(15) 


(16) 


(17) 


where 
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(18a) 


bt  =  c,  --  <£  %  -  «x|2)  A 

P  m  dp  '  dT  p' 

bp 

P  ,=I  op 


(ip)r  _(_^)r 

dp,.  ^dpN  Pi*N 


(18b) 

(18c) 


Enthalpy  as  function  of  temperature,  pressure,  and  mass  fraction 

The  third  relationship  expresses  enthalpy  as  a  function  of  pressure,  temperature,  and  mass 
fractions.  According  to  basic  thermodynamic  definitions,  we  have 


dh  -de  +  ~dp — ^ -dp 

p  p2 


(19) 


Substitution  of  Eq.  (8)  and  Eq.  (17)  into  Eq.  (19)  and  following  some  straightforward 
manipulations,  we  obtain 


A! -I 


dh  =  DTdT+  Dpdp  +  £  DY  dYi 


where 


Dr  =  C,  -  (£  Y,e,  -  e  - ?-) 

p  dT  p'  dp  •  tt  P 


P  P  dp  ‘  tt  p 


D,.=e, 

dp  i'=I  P 


dp 

' dPi 


U  h,plt,  )T,PjtN 


dp 

dpN 


(20) 


(21a) 


(21b) 

(21c) 


The  coefficient,  DT ,  is  equivalent  to  the  constant-pressure  heat  capacity,  Cp ,  of  a  fluid  mixture, 
according  to  its  definition.  Thus, 


C=DT=C--{dp 


dT 


dp 

dp 


(  N 


Jt.y, 


Yrti-e-z 

P) 


\  i-i 


(22) 


Relationship  between  constant-pressure  and  constant-volume  specific  heats 
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The  final  relationship  deals  with  constant-pressure  and  constant-volume  specific  heats. 
According  to  fundamental  thermodynamics  for  a  multi-component  mixture, 


s  =  s(T,p,Yi)  i  = 

Its  differential  form  can  be  written  as 


(23) 


ds=YdT  +  dP*%  <  JjT  W  dY‘ 


(24) 


A  similar  relationship  that  defines  entropy  as  a  function  of  temperature,  pressure,  and  mass 
fraction  can  be  easily  derived  as 

ds  =  ^-dT  + <^)r,n  dp  +  dr,  (25) 

Combination  of  Eqs.  (24)  and  (25)  and  rearrangement  of  the  result  lead  to 

Q-f)dT=(^)Trdp-(^)T,vdp  +  g[(|r>w,„  -<|rW  W  <26> 

Substitution  of  Eq.  (8)  into  Eq.  (26)  and  elimination  of  dp  give. 


r  Cp  C  05  dp  ds  ds  dp 

1t  t  W’-'W7* Vp 

.ds .  .05.  .dp.  ,  ... 

1=1  uij 


dr 


dp  '  dr/ 


(27) 


Since  temperature,  density,  and  mass  fraction  can  vary  independently,  the  coefficients  of  the 
above  differentials  must  vanish.  Thus, 


C  -C  -T(—)  (— ) 

'  v  {dpTY'dT)pY' 


(28) 


Application  of  the  well-known  Maxwell  relation  [21], 

A  =J_A  =__LA  /A 

{dPTY‘  PiKdT)pY‘  piKdT pYi rdpTYi’ 


(29) 


gives, 


C  =C  + 


11  „2 


dp_ 

\dTj 


V  ^ P  Jt,y> 


(30) 
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This  relationship  will  be  employed  to  calculate  the  eigenvalues  of  the  preconditioned  governing 
system. 


1.2.2  Preconditioning  Matrix 

The  precondition  matrix  r  in  Eq.  (2.1)  is  introduced  to  rescale  the  system  eigenvalues  in 
such  a  maimer  that  they  have  the  same  order  of  magnitude  to  ensure  uniform  convergence  at  all 
Mach  numbers.  The  specific  form  of  the  matrix  is  not  unique  [8],  and  is  selected  to  maximize  the 
numerical  efficiency  over  a  broad  range  of  flow  conditions.  Generally,  this  matrix  can  be 
constructed  by  means  of  either  an  eigenvalue  analysis  [3,8]  or  an  asymptotic  theory  [22,23]. 

Following  the  standard  approach  suggested  in  Refs.  [5],  we  examine  the  Jacobian  matrix 
T  =  dQ/dQ  in  Appendix  A.  A  common  term  (dp  / dp) T  y  ,  which  can  be  expressed  as 


PJt.y, 


r  f  \ 
S 


dp 


(2.31) 


is  identified  in  each  element  of  the  first  column  of  T.  It  is  this  term,  after  multiplication  of  the  time 
derivatives  of  dependent  variables  in  Eq.  (2.1),  that  dictates  the  propagation  speeds  of  acoustic 
waves  in  the  governing  system  [5].  The  preconditioning  matrix  is,  then,  established  by  simply 
replacing  this  term  with  a  scaling  factor,  0 . 


r= 


0 


0M 


0V 


4. 

AjU 

0  />  -¥■ 


0  0 


p  0 


P  dp 


©f, 


er*-, 


o  o 


0  0- 


\Yi 


Arjw-l 


AyU 


V 


04+(Z%-«-£K^)rj}  P“  Pv  P^r  ~~~~  pBr> 


p- 


\Yi 


\Y»-i 


Ay  U 

*N- 1 


Ay  V 


PBY  ~ 

'  *N- 1 


Ay  e, 

rN-\  1 


Ay  Y. 

rAT->  1 


P~ 


v  ) 


(32) 

where  the  two  sets  of  coefficients,  (Aj,  Ap,Ay)  and  ( BT ,  BY  )  are  defined  by  Eqs  (9)  and  (18), 
respectively.  The  scaling  factor,  0 ,  is  defined  as. 


0  =  [l/f2+(X-l)]/a2 


(33) 


The  preconditioning  factor,  e,  lies  between  0  and  1.  In  comparison  with  the 
preconditioning  matrices  defined  in  [4,5,13,15],  all  of  the  off-diagonal  terms  in  Eq.  (32)  are 
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retained.  The  preconditioning  matrix  becomes  the  Jacobian  matrix  T  =  dQ/dQ  in  the  limit  of 
£  — >  1 ,  and  the  original  governing  system  is  recovered  for  high-speed  flows  in  a  steady  state 
simulation.  For  an  ideal  gas,  the  preconditioning  matrix,  r ,  reduces  to  the  form  introduced  by 
Weiss  and  Smith  [6],  and  becomes  a  member  of  the  generalized  preconditioner  family  suggested  by 
Turkel  [3,8].  For  an  incompressible  fluid,  as  the  speed  of  sound,  a,  reaches  infinite,  the 
preconditioning  matrix  takes  a  form  similar  to  that  proposed  in  Ref.  [5].  As  will  be  shown  later, 
through  an  appropriate  selection  of  the  preconditioning  factor,  f ,  the  preconditioning  matrix 
developed  herein  can  effectively  circumvent  the  eigenvalue-disparity  problem  associated  with  low- 
Mach  number  flows  and  help  achieve  an  optimal  convergence  rate  in  the  entire  flowfield. 

According  to  Ref.  [5],  the  preconditioning  factor  for  an  inviscid  flow  is  specified  as 


f inv 


£min  M  <  £mm 


M 

1  M>  1 


(34) 


A  lower  bound,  fmin ,  typically  taken  as  10  5 ,  is  used  to  avoid  the  singularity  of  a  stagnation 
point.  In  regions  where  diffusion  plays  a  dominant  role  in  determining  the  flow  behavior,  it  is 
important  to  simultaneously  control  both  the  Courant-Friedriches-Lewy  (CFL)  and  von  Neumann 
(VNN)  numbers  to  achieve  efficient  convergence  [24].  Oefelein  and  Yang  [12]  suggested  that  the 
viscous  preconditioning  factor,  £vis ,  be  selected  as 


=  max 


u252x+a2  ’ 


v2Sy(Sy-l) 
v2 52  +  a2 


(35) 


where 


-  ,  v  v  ,1  CFL 

A  Pr  Sc/ u  VNN 


Sy  =  max(v', 


V_  v  1  CFL 
Pr’  Sc?  v  VNN 


(36) 


and  v ,  Pr ,  and  Sct  are  the  kinematic  viscosity,  cell  Prandtl  number,  and  cell  Schmidt  number  of 

species  i ,  respectively.  Equation  (36)  takes  into  account  the  effects  of  momentum,  energy,  and 
mass  diffusion  on  the  overall  convergence  rate. 

A  preconditioning  factor  that  can  optimally  control  both  convection  and  diffusion  processes 
is  determined  locally  as 

£  =  min[l,  max(£inv  ,£vis )]  (37) 


1.2.3  System  Eigenvalues 

The  basic  characteristics  of  the  algorithm  developed  herein,  specifically,  the  numerical 
convergence  and  stability  properties,  can  be  examined  by  studying  the  system  eigenvalues.  For 
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brevity,  only  those  eigenvalues  associated  with  the  inviscid  flux  vector  in  the  axial  direction,  T  'A 
with  A  =  dE /dQ  ,  are  discussed. 


Aa  =  diag(Ay  A^  u  u 


u) 


(38) 


where  \  and  A,  represent  the  rescaled  acoustic-wave  speeds  propagating  upstream  and 
downstream,  respectively. 


4.2  =\^n  +«)±V(^1  _M)2  -4(“^I1  -ea1)} 

The  variable  qn  represents  the  first  element  in  the  first  column  of  the  matrix,  T'A . 


(39) 


(K  t  g) 

( yc ,  T 

U2  p2K) 

! 

{/>  p'K) 

qu  —  u 


where  AT  and  A  ,  are  given  in  Eqs.  (9a)  and  (9c),  respectively.  We  then  obtain, 


(40) 


2  f  \2  ( 


AL_(dp^ 

ar 


dp) 

\  a/2  )t,y, 


(41) 


Combination  of  Eqs.  (41)  and  (30)  and  substitution  of  Eq.  (10)  into  the  resultant  equation  yield 

(42) 


T  4  C, 


C  =  C  '  +  ,  , 

'  p  K  cP 


v  a2 


The  above  equation  can  be  rearranged  to  become 

K  _  p 2  (r-Dcp 


a;  t  a i 


(43) 


Incorporation  of  Eq.  (43)  into  Eq.  (40)  leads  to 


qu  =U 


f rcr  t  a,2" 

/ 

3 

U2  p2K) 

/ 

{/>  p'K) 

=  eu 


(44) 


The  first  two  system  eigenvalues,  Ay2 ,  can  thus  be  obtained  by  substituting  Eq.  (44)  into  Eq.  (39). 

1, 


^  2  =  -[(1  +  £ )u  ±  yj(l-e)2u2  +4 £a2] 


(45) 


All  of  the  eigenvalues  in  the  pseudo-time  space  are  real,  and  have  signs  consistent  with  the 
directions  of  characteristic  wave  propagation.  The  present  scheme  not  only  preserves  the 
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hyperbolicity  of  the  system,  but  also  gives  rise  to  individual  eigenvalues  that  behave  in  a  manner 
representative  of  the  physical  reality  involved  [13].  The  same  conclusion  was  reached  in  [15]  using 
different  pseudo-time  variables  and  preconditioning  matrix. 

1.3  Numerical  Implementation 

A  dual-time-stepping  integration  technique  is  implemented  to  obtain  time-accurate  results 
[4,12,13,15],  The  solution  converged  in  pseudo-time  corresponds  to  a  time-accurate  solution  in 
physical-time.  One  major  advantage  of  this  technique  lies  in  the  fact  that  the  convergence  of  the 
iterative  process  is  dictated  by  the  well-behaved  eigenvalues  in  the  pseudo-time  space,  instead  of 
the  original  eigenvalues  that  may  become  disparate  in  certain  flow  regimes  (e.g.,  low  Mach-number 
flows)  [4], 

A  standard  fourth-order  Runge-Kutta  (RK4)  scheme  [17]  is  employed  to  perform  the  inner- 
loop  pseudo-time  integration  because  of  its  relatively  higher  temporal  accuracy  and  greater  stability 
margin  compared  to  many  commonly  used  explicit  schemes.  The  temporal  discretization  of  the 
real-time  derivative  term  is  obtained  using  a  second-order  backward  difference.  Spatial 
discretization  is  achieved  by  means  of  a  fourth-order  flux  differencing  scheme  [18].  Further 
improvement  is  acquired  by  adding  both  the  second-  and  fourth-order  artificial  dissipation  with  a 
total-variation-diminish  (TVD)  switch  [25,26]  to  ensure  numerical  stability  and  convergence.  The 
entire  approach  guarantees  that  the  second-order  artificial  dissipation  exerts  an  appropriate 
influence  in  regions  with  strong  gradients,  whereas  only  the  fourth-order  dissipation  exists  in 
smooth  regimes  to  achieve  numerical  stability.  The  details  regarding  implementation  of  artificial 
dissipation  in  a  preconditioning  scheme  can  be  found  in  [8]. 

Because  the  time  advancement  is  fully  explicit  in  pseudo-time  and  only  the  neighboring  data 
at  the  previous  time  step  is  required  when  evaluating  the  derivatives  of  the  convective  and  viscous 
terms,  the  present  scheme  is  suitable  for  parallel  implementation.  The  scheme  is  then  parallelized 
using  a  multiblock  domain  decomposition  technique  with  message  passing  interfaces  (MPI)  at  the 
domain  boundaries.  The  parallelization  methodology  is  robust  and  the  speedup  is  almost  linear. 

1.4  Stability  Analysis 

A  von  Neumann  analysis  is  performed  to  characterize  the  stability  and  convergence 
behavior  of  the  present  scheme.  For  simplicity,  only  the  inviscid  part  without  any  source  terms  is 
considered,  and  the  coefficient  matrices  are  assumed  to  be  locally  constant.  The  amplification 
matrix  is  obtained  by  the  Fourier  transformation  of  the  descretized  governing  equations  to  the 
wave-number  space.  The  partial  derivatives  of  the  thermodynamic  properties  in  the  corresponding 
numerical  Jacobian  matrices  and  the  fluid  p-v-T  behavior  are  evaluated  with  a  modified  Soave- 
Redlich-Kwong  equation  of  state  [27,28], 

Figure  la  shows  the  one-dimensional  stability  characteristics  for  an  ideal  gas  in  terms  of  the 
magnitude  of  the  largest  eigenvalue  of  the  amplification  matrix.  The  numerical  parameters  are 
CFL  =  0.7,  At/ At  =  1,  and  e[4>  =1/64,  where  the  latter  denotes  the  coefficient  of  the  fourth-order 

artificial  dissipation.  The  result  of  the  RK4  integration  combined  with  the  fourth-order  central 
differencing  (4CD)  scheme  is  also  included  for  comparison.  The  amplification  factor  of  the 
standard  RK4-4CD  algorithm  approaches  unity  in  the  bulk  of  the  wave-number  space,  whereas  the 
preconditioning  technique  stabilizes  numerical  calculations  over  the  entire  domain  at  all  Mach 
numbers. 
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The  effects  of  fluid  state  on  the  numerical  stability  behavior  are  investigated  by  considering 
nitrogen  fluid  at  a  near-critical  (p  =  40  atm  and  T  =  120  K  ),  a  transcritical  (p  =  90  atm  and 
T  =  \20K),  and  a  supercritical  ( p  =  90  atm  and  7  =  185  K)  fluid  state.  Table  1  lists  the 
thermodynamic  critical  properties  of  nitrogen.  The  compressibility  factors,  Z  ,  which  measure  the 
departure  from  the  ideal-gas  behavior,  are  0.2,  0.4,  and  0.8  at  those  three  different  fluid  states, 
respectively.  As  evidenced  Fig.  lb,  the  present  scheme  appears  to  be  quite  robust,  and  exhibits  an 
almost  identical  stability  behavior  independent  of  fluid  states.  The  same  result  is  obtained  for  all 
flow  Mach  numbers. 


Table  1  Critical  properties  of  nitrogen,  methane  and  oxygen. 


Pc  (atm) 

Tc  (K) 

vc  (L/mol) 

Nigrogen 

34 

126 

0.089 

Methane 

46 

190 

0.099 

Oxygen 

50 

154 

0.076 

Fig.  1  Amplification  factor  showing  the  one-dimensional  stability  characteristics  of  the  scheme 
at  CFL- 0.95,  At/ At  =  1 ,  e4  =  1/64 ,  a)  ideal  gas  at  M  =10~3,10“2,0.1,  and  1.25;  and  b) 
compressed  nitrogen  fluid  at  M  =  0.1  and  Z  =  0.2,  0.4,  0.8. 


Two-dimensional  stability  analyses  are  also  conducted  for  an  ideal  gas,  a  cryogenic  oxygen 
fluid  ( p  =  100  atm  and  T  =  100  K ),  and  a  supercritical  fluid  mixture  of  oxygen  and  methane 

( p  =  100  atm ,  T  =  200  K ,  and  xQ2  =  0.5 )  at  a  flow  Mach  number  of  10~3  and  a  flow  angle  of  45 
deg  (i.e.,  v/u  =  1).  The  compressibility  factors  for  the  latter  two  cases  are  0.34  and  0.48, 
respectively.  The  amplification  factors  shown  in  Fig.  2  indicate  that  the  scheme  exhibits  an 
identical  stability  behavior  in  a  two-dimensional  domain,  regardless  of  the  fluid  state.  The  result 
further  corroborates  the  numerical  uniformity  and  self-consistence  of  the  general  fluid 
thermodynamic  treatment  [15]  implemented  herein. 
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Z  =  1.0  Z  =  0.34  Z  =  0.48 


Fig.  2  Amplification  factors  showing  the  two-dimensional  stability  characteristics  of  the  scheme 
at  M  =  10“3 ,  CFL  =  0.95 ,  and  At/ At  =  1 ,  a)  ideal  gas,  Z  =  1 .0 ;  b)  compressed  cryogenic 
oxygen  at  100  K  and  100  atm,  Z  =  0.34 ;  c)  supercritical  oxygen/methane  mixture  at  200  K,  100 
atm,  and  oxygen  mole  fraction  of  0.5,  Z  =  0.48  . 


1.5  Sample  Calculations 

The  numerical  scheme  developed  in  the  preceding  sections  has  been  applied  to  study  a  wide 
variety  of  flow  problems  in  order  to  assess  its  accuracy,  efficiency,  and  robustness.  This  section 
presents  some  representative  results,  including  injection  of  cryogenic  fluids  and  mixing  and 
combustion  of  coaxial  oxygen/methane  fluid  jets  under  supercritical  conditions.  For  all  the 
demonstration  cases  considered  herein,  turbulence  closure  is  achieved  by  means  of  a  large-eddy- 
simulation  (LES)  technique,  in  which  large-scale  motions  are  calculated  explicitly  and  the  effects  of 
unresolved  small-scale  turbulence  are  modeled  either  analytically  or  empirically.  The  Favre- 
filtered  mass,  momentum,  energy,  and  species  conservation  equations  are  derived  by  filtering 
small-scale  dynamics  from  resolved  scales  over  a  well-defined  set  of  spatial  and  temporal  intervals. 
The  effects  of  subgrid-scale  (sgs)  motions  are  treated  using  the  model  proposed  by  Erlebacher  et  al. 
[29].  It  employs  a  Favre-averaged  generalization  of  the  Smagorinsky  eddy-viscosity  model  coupled 
with  a  gradient-diffusion  assumption  to  simulate  sgs  energy  and  species  transport  processes.  The 
Smagorinsky  coefficients  CR  (=0.01)  and  C,  (=0.007)  are  determined  empirically. 

Thermodynamic  properties,  such  as  enthalpy,  Gibbs  energy,  and  constant-pressure  specific  heat,  are 
obtained  directly  from  fundamental  thermodynamics  theories  and  a  modified  Soave-Redlich- 
Kwong  equation  of  state  [27,28],  Transport  properties,  such  as  viscosity  and  thermal  conductivity, 
are  evaluated  using  an  extended  corresponding-state  theory  [30,31]  along  with  the  32-term 
Benedict-Webb-Robin  equation  of  state  [32],  Mass  diffusivity  is  obtained  by  means  of  the 
Takahashi  method,  calibrated  for  high  pressure  conditions  [33].  The  implementation  and  validation 
of  the  property  evaluation  schemes  were  detailed  by  Yang  [1]  and  Meng  et  al.  [16]. 

1.5.1  Cryogenic  Nitrogen  Fluid  Jet  Dynamics 

The  first  case  deals  with  the  injection  of  cryogenic  fluids  into  supercritical  environments. 
Liquid  nitrogen  at  a  temperature  of  120  K  is  injected  through  a  circular  tube  with  a  diameter  of  254 
|im  into  a  supercritical  nitrogen  environment.  A  turbulent  pipe  flow  with  a  bulk  velocity  of  15  m/s 
is  assumed  at  the  injector  exit.  The  ambient  temperature  remains  at  300  K,  but  the  pressure  varies 
from  69  to  93  atm,  which  is  comparable  to  the  chamber  pressures  of  many  operational  rocket 
engines.  Two  different  flow  conditions  summarized  in  Table  2  are  considered,  simulating  the 
experiments  conducted  by  Chehroudi  and  colleagues  [34,35],  The  subscripts  «  and  inj  denote  the 
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injection  and  ambient  conditions,  respectively.  The  Reynolds  number  is  defined  as 
Re  =  p  u  D  .  /  u. .  .  . 

r  inj  inj  inj  f  r~uij 


Table  2  Simulation  conditions  for  injection  of  liquid  nitrogen. 


Poo 

(MPa) 

Tinj 

(K) 

Too 

(K) 

Pinj 

(kg/m3) 

Poo 

(kg/m3) 

Pinj/poo 

Uinj 

(m/s) 

Mjnj 

Re 

Case  1 

6.9 

jEEjl 

603 

mm 

44700 

Case  2 

9.3 

mSm 

626 

MM 

42300 

The  computational  domain  downstream  of  the  injector  measures  a  length  of  40  Dinj  and  a 
diameter  of  12  Dmj .  The  dimensions  are  sufficient  to  minimize  the  effect  of  the  far-field  boundary 

conditions  on  the  near-injector  flow  evolution.  A  three-dimensional  grid  consisting  of 
225x90x72  cells  is  employed.  The  mean  grid  size  falls  roughly  in  the  inertial  sub-range  of  the 
turbulent  kinetic  energy  (TKE)  spectrum,  estimated  using  the  Kolmogorov-Obukhow  theory.  The 
computational  domain  is  divided  into  54  blocks,  with  each  calculated  on  a  single  processor  of  a 
distributed-memory  parallel  computer.  The  physical  time  step  is  lxlO"3  ms  and  the  maximum 
CFL  number  for  the  inner-loop  pseudo-time  integration  is  0.7.  For  each  case,  the  calculation  is  first 
conducted  for  an  extended  period  until  the  flowfield  reaches  its  stationary  state.  The  time  stamp  is 
then  reset,  and  data  is  collected  for  more  than  12  flow-through  times  (i.e.,  15  ms)  to  obtain 
statistically  meaningful  turbulence  properties. 

Figure  3  shows  snapshots  of  the  density-gradient  fields  at  two  different  ambient  pressures  of 
6.9  and  9.3  MPa.  The  salient  features  of  supercritical  fluid  jets  are  well  captured.  The  fluid  state 
changes  continuously  from  the  injected  liquid  phase  to  the  warmer  ambient  gas  phase  with  a  series 
of  finger-  or  thread-like  entities  emerging  from  the  jet  surface  and  dissolving  gradually  in  the 
surrounding  gases.  Strong  anisotropy  of  turbulence  occurs  close  to  the  jet  interface,  where  large 


(b)  0  5  10  15  20 


Fig.  3  Effect  of  pressure  on  the  density-gradient  magnitude  field  (Too  =  300  K,  uinj  =  15  m/s,  T^j 
=  120  K,  Din,  =  254  Lim). 


Fig.  4  Radial  distributions  of  normalized  density  at  different  axial  locations  (T,*,=  300  K,  umj  = 

15  m/s,  T,nj  =  120  K,  Di„j  =  254  Jim). 

eddies  of  integral-length  scales  become  flattened,  and  the  radial  component  of  the  TKE  is 
transferred  to  its  axial  quantity  [36].  Compared  with  incompressible  turbulent  jets,  both  vortex  roll¬ 
up  and  pairing  are  delayed,  which  leads  to  a  longer  potential  core,  around  8 Dinj  for  Case  1.  The 
influence  of  the  density  stratification  decays  as  the  ambient  pressure  increases.  Thus,  the  location 
of  vortex  roll-up  shifts  upstream  from  x/ Dinj  ~  5  in  Case  1  to  i/ Dinj  ~  3  in  Case  2,  and  the  jet 

spreads  wider  and  the  length  of  the  potential  core  reduces  to  6-1  Dinj  in  the  latter  case. 

Figure  4  shows  the  radial  distributions  of  the  normalized  mean  density, 
p+  =(p- pm)l(pc  - pj),  at  different  axial  locations.  The  radial  coordinate  is  normalized  by  the 
full  width  of  the  radial  profile  measured  where  the  fluid  density  is  one  half  of  its  maximum  value 
(FWHM),  rU2,  at  the  axial  position  of  concern.  Three  distinct  flow  regimes,  similar  to 

incompressible  turbulent  jets,  are  identified.  The  potential  core  is  manifested  by  the  flat-hat 
distribution  near  the  injector.  The  density  profiles  merge  into  a  single  distribution  farther 
downstream  (x/ Dinj  >  30),  which  suggests  the  existence  of  a  fully  developed  self-similarity  region. 

A  transition  region  occurs  between  10  and  30 DinJ.  The  experimental  data  obtained  using  the 

Raman  scattering  technique  at  the  test  condition  of  Case  1  [34,35]  is  also  plotted  for  comparison. 
Good  agreement  between  calculations  and  measurements  is  achieved  with  a  maximum  deviation  of 
8%. 


1.5.2  Supercritical  Mixing  and  Combustion  of  Oxygen  and  Methane 

The  second  case  treats  the  flow  and  flame  dynamics  of  a  shear-coaxial  injector  operating  at 
supercritical  conditions,  as  shown  schematically  in  Fig.  5.  Co-flowing  methane  and  oxygen  streams 
are  injected  through  two  concentric  tubes.  The  inner  diameters  of  the  LOX  post  and  the  methane 
annulus  are  3.42  and  5.18  mm,  respectively.  Fully  developed  turbulent  pipe  flows  are  assumed  at 
the  injector  exit.  In  the  cold-flow  simulations,  supercritical  oxygen  and  methane  are  injected  at 
temperatures  of  200  and  300  K,  respectively.  The  bulk  velocities  of  the  two  streams  are  10  and  30 
m/s,  respectively. 

The  computational  domain  extends  AOS  downstream  of  the  injector  exit  with  a  radius  of 
12 S,  where  8  is  the  thickness  of  the  LOX  post.  A  quasi-axisymmetric  simulation  involving 
360x200  grid  pints  was  conducted.  The  mean  grid  size  of  10  pm  is  sufficient  to  resolve  the 
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Fig.5  Schematic  diagram  of  shear  coaxial  injection  of  oxygen  and  methane. 

inertial  sub-range  of  the  TKE  spectrum  estimated  based  on  the  inlet  Reynolds  number  of  the 
oxygen  stream.  The  physical  time  step  is  0.5  xlCT3  ms  and  the  maximum  CFL  number  for  the 
inner-loop  pseudo-time  integration  is  0.7.  A  more  detailed  description  of  the  simulation  conditions 
is  given  in  [37]. 


Figure  6  presents  close-up  views  of  the  vorticity,  temperature,  oxygen  mass-fraction,  and 
compressibility-factor  fields  near  the  injector  without  chemical  reactions.  Three  shear  layers  are 
clearly  observed:  one  emerging  from  the  inner  rim  of  the  LOX  post;  and  two  from  the  inner  and 
outer  edges  of  the  methane  annulus.  A  series  of  large  scale  vortices  shed  from  the  outer  rim  of  the 
LOX  post.  As  those  vortices  develop,  the  two  shear  layers  separated  by  the  LOX  post  merge. 
Although  those  energetic  eddies  concentrate  on  the  light  fluid  side,  they  entrain  the  denser  oxygen 
much  deep  into  the  methane  stream  and  greatly  facilitate  the  mixing  process. 

The  situation  with  chemical  reactions  is  also  considered.  Liquid  oxygen  (LOX)  and 
methane  are  injected  at  temperatures  of  122  and  300  K,  respectively.  The  bulk  velocities  of  the  two 
streams  are  13  and  75  m/s,  respectively.  The  mixture  and  momentum-flux  ratios  are  3  and  2.5, 
respectively.  The  injection  conditions  are  typical  of  operational  liquid-propellant  rocket  engines. 
The  combustion  chamber  is  preconditioned  with  a  mixture  of  COz  and  H20  at  the  stoichiometric 
ratio  and  1800  K.  The  ambient  pressure  is  fixed  at  100  atm.  A  one-step  global  chemical  kinetics 
model  for  methane  and  oxygen  is  employed  [38]. 


Figure  7  shows  snapshots  of  vorticity,  temperature,  methane  mass-fraction,  and 
compressibility-factor  fields.  A  diffusion-dominated  flame  emanates  immediately  from  the  LOX 
post  and  propagates  downstream  along  the  surface  of  the  LOX  stream.  A  wake  region,  which 
consists  of  hot  combustion  products,  effectively  separates  the  methane  and  LOX  streams.  Similar 
to  the  non-reacting  flow  case,  the  near-field  flow  dynamics  are  characterized  by  the  evolution  of  the 
three  mixing  layers  originating  from  the  inner  and  outer  edges  of  the  methane  annulus  and  the  inner 
rim  of  the  LOX  post.  The  evolution  of  the  inner  mixing  layer  of  the  methane  stream,  however,  is 
slightly  inhibited  by  the  expansion  of  the  combustion  products  in  the  flame  zone.  Because  of  the 
strong  density  stratification  between  the  oxygen  stream  and  flame,  the  large-scale  vortices  emerging 
from  the  outer  rim  of  the  LOX  post  evolve  in  a  manner  analogous  to  that  produced  by  a  backward¬ 
facing  step  and  mainly  reside  on  the  lighter  fluid  side.  Consequently,  the  denser  oxygen  stream  is 
less  influenced  by  those  vortices. 
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Fig.  6  Snapshots  of  vorticity,  temperature,  oxygen  mass-fraction,  and  compressibility-factor 
fields  of  shear  coaxial  injection  of  oxygen  and  methane  (p,»  =  100  atm,  T02  =  200  K,  U02  =  10 
m/s,  Tch4=  300  K,  uch4=  30  m/s). 


To  assess  the  numerical  performance  of  the  present  scheme,  calculations  are  also  conducted 
using  the  preconditioning  scheme  developed  in  [15]  for  both  the  non-reacting  and  reacting  flows. 
The  same  flow  solver  described  in  Section  3  is  employed,  but  with  the  preconditioning  formulation 
based  on  the  pressure-enthalpy  type  of  primitive  variables  and  a  different  preconditioning  matrix  in 
[15].  The  use  of  enthalpy  as  one  of  the  dependent  variables  leads  to  an  iteration  solution  of  the  real- 
fluid  equation  of  state  by  means  of  the  Newton-Raphson  method  with  a  convergence  criterion  of 
|A/?|/p  <  lO^and  |AT|/r  <1(T*.  A  relaxation  of  this  criterion  may  result  in  an  overflow  of  the 

calculation  due  to  strong  property  variations  in  the  flowfield.  Figure  8  compares  the  convergence 
histories  of  computations  using  the  schemes  developed  in  the  present  study  and  in  [15].  Because 
the  same  temporal  and  spatial  discretization  techniques  are  employed,  the  two  schemes  exhibit  an 
almost  identical  behavior  of  convergence  in  terms  of  the  number  of  pseudo-time  iterations.  The 
present  scheme,  however,  is  much  more  efficient  because  it  avoids  laborious  iterations  in 
determining  temperature  from  enthalpy.  Table  3  lists  the  CPU  times  per  pseudo-time  iteration,  t  , 

and  the  times  for  solving  the  equation  of  state,  teos .  All  the  computations  were  performed  on  a 

Pentium  IV  2.4  GHz  processor.  Less  than  5%  of  the  total  CPU  usage  is  expended  for  solving  the 
equation  of  state  with  the  present  scheme,  as  opposed  to  more  that  50%  in  the  scheme  in  [15].  The 
overall  computational  time  is  also  reduced  by  half.  Another  major  benefit  of  the  present  scheme  is 
that  the  computational  burden  for  each  spatial  grid  cell  is  the  same  because  no  iterations  are 
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required  for  solving  the  equation  of  state.  The  load  balance  among  numerical  blocks  is  much  easier 
to  achieve  on  a  distributed  computing  facility. 


Fig.  7  Snapshots  of  vorticity,  temperature,  methane  mass-fraction,  and  compressibility-factor 
fields  near  the  injector  faceplate  (p«,  =  100  atm,  T02  =  122  K,  U02  =  13  m/s,  Tch4  =  300  K,  uch4  = 
75  m/s). 


Fig.  8  Convergence  histories  in  the  pseudo-time  domain  for  shear  coaxial  injection  and 
combustion  of  oxygen  and  methane. 
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Table  3  CPU  time  per  iteration  in 

pseudo-time. 

non-reacting  flow  case 

reacting  flow 

case 

*sp  00  Kos  00  t'oJhp 

hP  oo  teos  (s) 

t  It 

eos  sp 

Present  scheme 

52  1.8  3% 

87  3.8 

4% 

Meng  &  Yang  [15] 

120  70  58% 

158  86 

54% 

1.6  Conclusions 

A  general  treatment  of  real-fluid  thermodynamics  within  the  framework  of  a 
preconditioning  scheme  has  been  established.  All  of  the  numerical  properties  are  derived  directly 
from  fundamental  thermodynamics  theories  based  on  the  concepts  of  partial-mass  and  partial- 
density  properties.  The  algorithm  is  self-consistent  and  robust,  and  can  accommodate  any  equation 
of  state  for  fluid  mixtures.  The  scheme  employs  temperature  instead  of  enthalpy  as  the  primary 
dependent  variable  in  the  preconditioned  energy  equation.  As  a  consequence,  the  laboriously 
iterative  solution  of  the  equation  of  state  can  be  avoided.  The  computational  burden  is  uniform 
throughout  the  entire  domain.  A  series  of  sample  calculations,  including  the  injection  and 
combustion  of  cryogenic  fluids  under  supercritical  conditions,  were  conducted  to  assess  the 
effectiveness  of  the  scheme  at  various  fluid  states.  In  addition,  a  numerical  stability  analysis  was 
performed  to  characterize  the  algorithm  behavior  over  a  broad  range  of  flow  conditions. 
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Appendix  A 


The  Jacobian  matrices,  T  =  dQ/dQ ,  A  =  dEI dQ ,  B  =  dF  I  dQ ,  and  D  =  dS  Id Q  are  derived  as 
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Task  2 

Cryogenic  Propellant  Injection  and  Mixing  at  Supercritical  Conditions 


This  work  treats  cryogenic  propellant  injection  and  mixing  at  near-critical  and  super-critical 
conditions  under  conditions  with  and  without  transverse  acoustic  excitations,  simulating  the 
experiments  conducted  by  Chehroudi  and  Talley  at  the  Air  Force  Research  Laboratory.  The 
purposes  are  to  achieve  an  improved  understanding  of  high-pressure  cryogenic  fluid  jet  dynamics, 
and  to  further  validate  the  code  outlined  in  the  subtask  1.  The  analysis  accommodates  full 
conservation  laws  and  real-fluid  thermodynamics  and  transport  over  the  entire  temperature  and 
pressure  regimes  of  concern.  Turbulence  closure  is  achieved  using  a  large-eddy-simulation  (LES) 
technique,  in  which  large  energy-carrying  structures  are  computed  explicitly  and  the  effects  of 
small-scale  motions  on  resolved  scales  are  modeled.  A  Smagorinsky  model  extended  to 
compressible  flows  is  used  to  treat  small-scale  motions.  Thermodynamic  properties,  including 
enthalpy,  internal  energy,  and  heat  capacity,  are  directly  calculated  by  means  of  fundamental 
thermodynamic  theories  and  a  modified  Soave-Redlich-Kwong  (SRK)  equation.  Transport 
properties,  such  as  viscosity  and  thermal  conductivity,  are  estimated  with  an  extended 
corresponding-state  principle.  Mass  diffusivity  is  obtained  by  the  Takahashi  method  calibrated  for 
high-pressure  conditions.  Numerical  implementation  includes  a  preconditioning  scheme  in 
conjunction  with  a  dual  time-stepping  integration  algorithm.  Further  efficiency  is  achieved  by 
utilizing  a  parallel  computational  scheme  that  involves  the  message-passing-interface  (MPI)  library 
and  a  multi-block  treatment.  Results  from  this  task  have  led  to  one  archival  journal  publication  as 
follows. 


“A  Numerical  Study  of  Cryogenic  Fluid  Injection  and  Mixing  under  Supercritical 
Conditions,”  by  N.  Zong,  H.  Meng,  S.Y.  Hsieh,  and  V.  Yang,  Physics  of  Fluids,  Vol.  16, 
2004,  pp.  4248-4261. 

2.1  Introduction 

Injection  of  liquid  fluid  initially  at  a  subcritical  temperature  into  an  environment  in  which 
the  temperature  and  pressure  exceed  the  thermodynamic  critical  point  of  the  injected  fluid  is  an 
important  phenomenon  in  many  high-pressure  combustion  devices,  including  diesel,  gas-turbine, 
and  liquid-propellant  rocket  engines.  Under  these  conditions,  fluid  jets  exhibit  many  characteristics 
distinct  from  their  counterparts  at  low  pressures.  First,  because  of  the  diminishment  of  surface 
tension  and  enthalpy  of  vaporization,  the  sharp  distinction  between  the  liquid  and  gas  phases 
vanishes.  The  fluid  properties  and  their  spatial  gradients  vary  continuously  throughout  the  entire 
field.  Second,  thermodynamic  and  transport  anomalies  may  occur  during  the  temperature  transition 
across  the  critical  value,  especially  when  the  pressure  approaches  the  critical  point,  a  phenomenon 
commonly  referred  to  as  near-critical  enhancement  [1].  As  a  result,  compressibility  effects  (i.e., 
volumetric  changes  induced  by  pressure  variations)  and  variable  inertial  effects  (i.e.,  volumetric 
changes  induced  by  heat  addition  and/or  variable  composition)  play  an  important  role  in  dictating 
the  flow  evolution.  Third,  the  characteristic  times  of  the  flow  motions  around  the  jet  boundary  have 
the  same  order  of  magnitude  at  supercritical  conditions.  The  resultant  coupling  dynamics  becomes 
transient  in  both  the  jet  interior  and  the  surrounding  fluid,  and  involves  an  array  of  physiochemical 
processes  with  widely  disparate  time  and  length  scales.  Fourth,  the  flow  Reynolds  number 
increases  almost  linearly  with  pressure.  For  oxygen  and  hydrogen,  an  increase  in  pressure  from  1  to 
100  atmospheres  results  in  a  reduction  of  the  corresponding  kinematic  viscosity  by  a  factor  of  two 
orders  of  magnitude.  Based  on  the  Kolmogorov  universal  equilibrium  theory,  the  Kolmogorov  and 
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variable-density,  single-phase  turbulent  submerged  jet  at  both  subcritical  and  supercritical 
pressures,  when  the  ambient  temperature  remained  supercritical. 

Motivated  by  the  development  of  high-pressure  cryogenic-propellant  rocket  engines, 
extensive  experimental  studies  have  recently  been  conducted  to  provide  direct  insight  into 
supercritical  fluid  injection  and  mixing.  The  work  included  injection  of  liquid  nitrogen  and  co¬ 
injection  of  liquid  nitrogen  and  gaseous  helium  into  gaseous  nitrogen  environments  over  a  wide 
range  of  pressure  [8-10].  Cryogenic  propellants  such  as  liquid  oxygen  and  hydrogen  under  both 
cold-flow  and  hot-fire  test  conditions  were  also  considered.  Results  from  the  shadowgraph  images 
confirmed  the  findings  by  Newman  and  Brzustowski  [7].  Drastic  changes  in  jet  surface  phenomena 
took  place  across  the  critical  pressure.  Ligaments  and  droplets  formed  at  subcritical  pressures,  but 
disappeared  at  supercritical  conditions  due  to  the  prevalence  of  turbulent  motions  and  vanishment 
of  surface  tension.  The  jet  surface  topology  bears  a  strong  resemblance  to  its  gaseous  counterpart, 
with  the  spatial  growth  rate  following  that  of  an  imcompressible  variable-density  gaseous  jet  [10]. 
The  spontaneous  Raman  Scattering  technique  was  employed  to  measure  the  density  field  [11,12], 
In  general,  the  normalized  density  profiles  indicated  a  tendency  towards  the  self-similarity  solution 
observed  for  classical  constant-  and  variable-density  single-phase  fluid  jets  at  low  pressures.  The 
effects  of  acoustic  waves  on  cryogenic  nitrogen  jet  flow  development  were  recently  studied  by 
Chehroudi  and  Talleyf  13]  The  influence  was  substantial  at  subcritical  conditions,  but  became 
unnoticeable  at  supercritical  conditions.  The  phenomenon  may  be  attributed  to  the  formation  of 
high-frequency  vortices  in  the  supercritical  jet  flow. 

In  parallel  to  experimental  studies,  attempts  were  made  both  theoretically  and  numerically 
to  explore  the  underlying  mechanisms  of  high-pressure  fluid  injection  and  mixing.  Oefelein  and 
Yang  [2]modeled  two-dimensional  mixing  and  combustion  of  oxygen  and  hydrogen  streams  at 
supercritical  conditions  by  means  of  a  large-eddy-simulation  technique.  The  formulation 
accommodated  real-fluid  thermodynamics  and  transport  phenomena.  All  the  thermophysical 
properties  were  evaluated  directly  from  fundamental  thermodynamics  theories  over  the  entire 
regime  of  fluid  states  of  concern.  Furthermore,  a  unified  treatment  of  numerical  algorithms  based 
on  general  fluid  thermodynamics  was  established  to  improve  computational  accuracy  and 
efficiency[3,14]  Bellan  and  colleagues  [15,16,17]  treated  the  temporal  evolution  of 
heptane/nitrogen  and  oxygen/hydrogen  mixing  layers  at  supercritical  conditions.  Several  important 
characteristics  of  high-pressure  transitional  mixing  processes  were  identified.  In  particular,  as  a 
result  of  density  stratification,  the  mixing  layer  is  considerably  more  stable  than  a  corresponding 
gaseous  mixing  layer.  Energy  dissipation  due  to  both  the  species-flux  and  heat-flux  effects  prevails 
during  the  evolution  of  the  mixing  layer,  whereas  the  viscous  effect  appears  minimal. 

The  purpose  of  the  present  work  is  to  conduct  a  comprehensive  analysis  of  cryogenic  fluid 
injection  and  mixing  under  supercritical  conditions.  The  theoretical  formulation  accommodates  full 
conservation  laws  and  real-fluid  thermodynamics  and  transport  phenomena.  Turbulence  closure  is 
achieved  using  a  large-eddy-simulation  technique.  As  a  specific  example,  the  dynamics  of  a 
cryogenic  nitrogen  jet  in  a  supercritical  nitrogen  environment  is  studied  over  a  wide  range  of 
pressure.  The  work  focuses  on  the  near-field  jet  dynamics  and  its  subsequent  evolution. 

2.2  Theoretical  Formulation 

A  unified  theoretical  framework  capable  of  treating  fluid  flows,  trans-critical  property 
variations  and  real-fluid  thermodynamics  is  developed.  Turbulence  closure  is  achieved  using  a 
large-eddy-simulation  (LES)  technique,  in  which  large-scale  motions  are  calculated  explicitly  and 
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the  effects  of  unresolved  small-scale  turbulence  are  modeled  either  analytically  or  empirically. 
Thermodynamic  properties  such  as  enthalpy,  Gibbs  energy,  and  constant-pressure  specific  heat  are 
obtained  directly  from  fundamental  thermodynamics  theories  [3,14]  Transport  properties  are 
evaluated  using  an  extended  corresponding-state  theory  along  with  the  32-term  Benedict-Webb- 
Robin  equation  of  state  [3]. 


2.2.1  Governing  Equations 

The  formulation  is  based  on  the  Favre-filtered  conservation  equations  of  mass,  momentum 
and  energy.  These  equations  are  derived  by  filtering  the  small-scale  dynamics  from  the  resolved 
scales  over  a  well-defined  set  of  spatial  and  temporal  intervals.  They  can  be  conveniently 
expressed  in  the  following  Cartesian  tensor  form: 
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where  overbars  and  tildes  denote  resolved-scale  and  Favre-averaged  resolved-scale  variables, 
respectively,  p,  u„  p,  E,  and  q,  represent  the  density,  velocity  components,  pressure,  specific  total 
energy,  viscous  stress  tensor,  heat  flux,  respectively.  A  detailed  derivation  of  the  filtered  equations 
is  presented  by  Oefelein  and  Yang  [2].  The  unclosed  subgrid-scale  (sgs)  terms  in  Eqs.  (1)  -  (3), 
including  the  stresses  T-f” ,  energy  fluxes  H-gs  and  viscous  work  c'f ,  are  defined  as 
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The  unclosed  sgs  stress  term  is  modeled  using  the  compressible  flow  version  of  the 
Smagorinsky  model  suggested  by  Erlebacher  et  al.  [18] 

r;r  =-2r,p(5,-SuS,/3)  +  lfk-’S,  (7) 


where 


The  dimensionless  constants  Cr  (=0.01)  and  C/  (=0.007)  are  determined  empirically, 
ksgs  =ts£  !2p  =  (pukuk  /p-ukuk)/2  and  Sy  =(dui/dxj+duj/dxi)/ 2.  The  symbol  A  denotes 

the  filter  width.  The  subgrid  energy  flux  term  H.gs  is  modeled  as 
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A  standard  value  of  0.7  is  used  for  the  turbulent  Prandtl  number.  The  unclosed  subgrid  viscous- 
work  term  rx/f  is  assumed  to  be  small  and  is  neglected  in  the  present  study. 


2.2.2  Thermodynamic  and  Transport  Properties 

Owing  to  the  continuous  variations  of  fluid  properties  throughout  the  jet  and  surroundings 
under  supercritical  conditions,  conventional  treatment  of  fluid  jets  at  low  pressures  in  which  the 
liquid  and  gas  phases  are  solved  separately  and  then  matched  at  the  interfaces  often  leads  to 
erroneous  results.  The  problem  becomes  even  more  exacerbated  when  the  fluid  approaches  its 
critical  state,  around  which  fluid  properties  exhibit  anomalous  sensitivities  with  respect  to  local 
temperature  and  pressure  variations.  Thus,  a  prerequisite  of  any  realistic  treatment  of  supercritical 
fluid  behavior  lies  in  the  establishment  of  a  unified  property  evaluation  scheme  valid  over  the  entire 
thermodynamic  regime  [3,14]. 

The  formulation  established  here  can  accommodate  any  equation  of  state.  In  the  present 
work,  a  modified  Soave-Redlick-Kwong  (SRK)  equation  of  state  is  adopted  because  of  its  wide 
range  of  validity  and  ease  of  implementation  [19].  It  takes  the  following  form, 
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where  Ru  is  the  universal  gas  constant.  The  parameters  a  and  b  account  for  the  effects  of  attractive 

and  repulsive  forces  between  molecules,  respectively.  The  third  parameter  a  includes  the  critical 
compressibility  factor  and  acentric  (size-shape)  interactions  between  molecules. 


Thermodynamic  properties,  such  as  enthalpy,  internal  energy,  and  constant-  pressure 
specific  heat,  can  be  expressed  as  the  sums  of  ideal-gas  properties  at  the  same  temperature  and 
departure  functions  taking  into  account  the  dense-fluid  correction  [3].  Thus, 
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The  subscript  0  refers  to  the  ideal  state  at  low  pressure.  The  departure  functions  on  the 
right-hand  sides  of  equations  (12)-(14)  can  be  determined  from  the  equation  of  state. 

Transport  properties  including  viscosity  and  thermal  conductivity  are  estimated  by  means  of 
the  corresponding  state  principles  along  with  the  use  of  the  32-term  Benedict-Webb-Rubin(BWR) 
equation  of  state  [20,21].  Figure  2  shows  the  effects  of  temperature  on  the  compressibility  factor, 
specific  heat,  viscosity,  and  thermal  conductivity  of  nitrogen  in  the  temperature  range  of 
40  - 1000  K  at  five  different  pressures  of  1,  10,  33.4,  100,  200  atmospheres.  These  figures  illustrate 
key  trends  over  the  range  of  pressure  and  temperature  of  interest.  At  750  K  and  above,  nitrogen 
exhibits  ideal  gas  behavior  and  the  pressure  effect  is  negligible.  As  the  temperature  decreases 
below  750  K,  however,  significant  non-idealities  are  introduced  with  strong  pressure  dependence. 
The  implementation  and  validation  of  the  property  evaluation  schemes  were  discussed  in  Ref.  [3]. 
It  should  be  noted  that  the  present  scheme  does  not  consider  the  critical  enhancement  in  the  close 
vicinity  of  the  critical  point,  a  region  known  as  the  critical  region  proper  where  a  scaled  equation  of 
state  must  be  used  [22].  Thus,  the  thermophysical  properties  shown  in  Fig.  2  do  not  diverge  at  the 
critical  point  as  predicted  by  classical  theories.  Since  the  simulation  conditions  chosen  herein 
deviate  considerably  from  the  critical  point,  the  neglect  of  critical  divergence  exerts  no  influence  on 
the  solution  accuracy. 


■*  40  250  500  750  40  250  500  750 

Temperature,  K  Temperature,  K 

Fig.2  Thermo-physical  properties  of  nitrogen  as  functions  of  pressure  and  temperature. 
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2.3  Numerical  Method 

The  present  work  bears  two  severe  numerical  challenges.  First,  thermodynamic  non¬ 
idealities  and  transport  anomalies  take  place  as  the  fluid  transits  through  the  transcritical  regime. 
Thus,  treating  these  phenomena  in  a  manner  consistent  with  the  intrinsic  characteristics  of  a 
numerical  algorithm  presents  a  major  obstacle.  Second,  the  rapid  variation  of  fluid  state  and  wide 
disparities  in  the  characteristic  time  and  length  scales  pose  the  well-known  stiffness  problem.  A 
unified  treatment  of  general  fluid  thermodynamics  is  developed  to  circumvent  those  difficulties. 
The  method  is  based  on  the  concepts  of  partial-mass  and  partial-density  properties  [14],  and  valid 
for  the  entire  pressure  and  temperature  regimes  of  concern.  The  resultant  routine  is  incorporated 
into  a  preconditioning  scheme  along  with  a  dual  time-step  integration  algorithm  to  handle 
compressible  fluid  flows  at  all  speeds[23,24].  Since  all  of  the  numerical  relations,  including  the 
Jacobian  matrices  and  eigenvalues,  are  derived  directly  from  fundamental  thermodynamics  theories, 
the  algorithm  is  self-consistent  and  robust.  A  detailed  description  of  the  overall  approach  was  given 
by  Meng  and  Yang  [14]. 

The  theoretical  formulation  outlined  above  is  solved  by  means  of  a  density-based,  finite 
volume  methodology.  The  spatial  discretization  employs  a  fourth-order,  central-difference  scheme 
in  generalized  coordinates.  A  fourth-order  scalar  dissipation  with  a  total-variation-diminishing 
switch  developed  by  Swanson  and  Turkel  [25]  is  implemented  to  ensure  computational  stability  and 
to  prevent  numerical  oscillations  in  regions  with  steep  gradients.  Temporal  discretization  is 
obtained  using  a  second-order  backward  difference,  and  the  inner-loop  pseudo-time  term  is 
integrated  with  a  three-step  third-order  Runge-Kutta  scheme.  A  multiblock  domain  decomposition 
technique,  along  with  static  load  balance,  is  employed  to  facilitate  the  implementation  of  parallel 
computation  with  message  passing  interfaces  at  the  domain  boundaries. 

2.4  Computational  Domain  and  Grid  System 

The  physical  model  of  concern  is  shown  schematically  in  Fig.  1,  simulating  the  experiments 
described  in  Refs.  [10]  and  [12].  Cryogenic  nitrogen  fluid  is  injected  into  supercritical  gaseous 
nitrogen  through  a  circular  duct  with  an  inner  diameter  of  0.254  mm.  Because  of  the  enormous 
computational  effort  required  for  calculating  the  flowfield  in  the  entire  three-dimensional  regime, 
only  a  cylindrical  sector  with  periodic  boundary  conditions  specified  in  the  azimuthal  direction  is 
treated  herein.  The  analysis,  in  spite  of  the  lack  of  vortex-stretching  mechanism,  has  been  shown  to 
be  able  to  capture  the  salient  features  of  supercritical  fluid  jet  evolution  [26]. 

The  computational  domain  downstream  of  the  injector  measures  a  length  of  40  Dinj  and  a 
radius  of  6  Dmj .  The  dimensions  are  sufficient  to  minimize  the  effect  of  the  far-field  boundary 

conditions  on  the  near-injector  flow  evolution.  The  entire  grid  system  consists  225x90  points 
along  the  axial  and  radial  directions,  respectively.  The  grids  are  clustered  in  the  shear  layer  and 
near  the  injector  to  resolve  rapid  property  variations  in  those  regions.  The  mean  grid  size  in  the 
near  field  (0  <  x / Dinj  <  20)  falls  in  the  inertial  sub-range  of  the  turbulent  kinetic  energy  spectrum, 

estimated  using  the  Kolmogorov-Obukhow  theory.  The  smallest  grid  width  is  2  /mi .  The 
computational  domain  is  divided  into  45  blocks,  with  each  calculated  on  a  single  processor  of  a 
distributed-memory  parallel  computer.  The  physical  time  step  is  lxl0~3  ms  and  the  maximum 
CFL  number  for  the  inner-loop  pseudo-time  integration  is  0.7.  For  each  case,  simulation  was 
conducted  for  12  flow-through  times  (i.e.,  15  ms)  to  obtain  statistically  meaningful  data. 
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A  grid  independence  study  was  performed  as  part  of  the  validation  procedure,  in  which  the 
same  numerical  code,  configuration,  and  flow  condition  (i.e.  Case  1  in  Table  1)  were  considered 
with  two  different  grid  resolutions:  a  fine  (270x120)  and  a  coarse  (225x90)  mesh.  Figure  3 
shows  the  radial  distributions  of  the  mean  axial  velocity,  turbulent  kinetic  energy  (TKE), 
compressibility  factor  Z,  and  viscosity  (normalized  by  the  value  of  the  ambient  gas)  at  different 
axial  locations.  Results  from  the  two  different  grid  systems  agree  well  with  each  other,  except  for 
the  small  deviation  of  TKE. 


90x225  .  120x270 


Fig.3  Effect  of  grid  resolution  on  radial  distributions  of  mean  axial  velocity,  turbulent  kinetic 
energy,  compressibility  factor,  and  viscosity  at  different  axial  locations  (p=*>  =  4.3  MPa,  Too  = 
300  K,  Uinj  =  15  m/s,  Ti„j  =  120  K,  Dinj  =  254  (im) 


2.5  Boundary  Condition 

At  the  injector  exit,  a  fully  developed  turbulent  pipe  flow  is  assumed.  The  mean  velocity 
follows  the  one-seventh-power  law  and  the  temperature  is  specified  with  a  top-hat  profile.  The 
pressure  is  determined  using  a  one-dimensional  approximation  to  the  momentum  equation  in  the 
axial  direction.  Turbulence  is  provided  by  superimposing  broad -band  white  noise  onto  the  mean 
velocity  profile.  The  disturbances  are  generated  by  a  Gaussian  random-number  generator  with  an 
intensity  of  12%  of  the  mean  quantity.  At  the  downstream  boundary,  extrapolation  of  primitive 
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variables  from  the  interior  may  cause  undesired  reflection  of  waves  propagating  into  the 
computational  domain.  Thus,  the  non-reflecting  boundary  conditions  proposed  by  Poinsot  and  Lele 
[27]  are  applied,  along  with  the  specification  of  a  reference  pressure.  Because  the  jet  flow  is 
directly  exhausted  to  an  ambient  condition,  the  surrounding  fluid  may  be  entrained  into  the 
computational  domain.  At  the  radial  boundary,  the  pressure,  temperature,  and  axial  velocity  are 
specified.  The  conservation  law  of  mass  is  employed  to  determine  the  radial  velocity.  Finally,  the 
non-slip  adiabatic  conditions  are  enforced  along  the  solid  walls. 

2.6  Results  and  Discussions 

The  theoretical  model  and  numerical  scheme  established  in  the  preceding  sections  were 
implemented  to  study  the  injection  and  mixing  of  cryogenic  fluid  under  supercritical  conditions.  As 
a  specific  example,  liquid  nitrogen  at  a  temperature  of  120  K  is  injected  through  a  circular  tube  with 
a  diameter  of  254  pm  into  a  supercritical  nitrogen  environment.  A  turbulent  pipe  flow  with  a  bulk 
velocity  of  15  m/s  is  assumed  at  the  injector  exit.  The  ambient  temperature  remains  at  300  K,  but 
the  pressure  varies  from  42  to  93  atm,  comparable  to  the  chamber  pressures  of  many  operational 
rocket  engines.  For  reference,  the  critical  temperature  and  pressure  of  nitrogen  are  126  K  and  34 
atm,  respectively.  Three  different  flow  conditions  summarized  in  Table  1  are  considered, 
simulating  the  experiments  conducted  by  Chehroudi  and  Talley  [10],where  the  subscripts  °°  and  inj 
denote  the  injection  and  ambient  conditions,  respectively.  The  Reynolds  number  is  defined  as 

Re  =  P,njUinjDlr,j  /Min]  • 


Table  1  Simulation  Conditions 


Case  1  Case  2  Case  3 


Poo  (MPa) 

4.3 

6.9 

9.3 

Too  (K) 

300 

300 

300 

Poo 

46 

77 

103 

Tinj  (K) 

120 

120 

120 

Pinj 

563 

603 

626 

Uini  (m/s) 

15 

15 

15 

Pinj/ po<= 

12.24 

7.83 

6.07 

Re 

48500 

44700 

42300 

Figure  4  shows  the  variations  of  nitrogen  density  and  constant-pressure  specific  heat  as 
functions  of  temperature  at  four  different  pressures.  Two  observations  are  noted  here.  First,  the 
density  decreases  sharply  near  the  critical  point  as  the  temperature  increases.  The  effect  of  density 
stratification  between  the  jet  and  the  ambient  fluid  becomes  much  more  substantial  for 
P"  =  4.3  MPa  compared  with  the  other  two  cases.  Second,  the  temperature  sensitivity  of  the 
specific  heat  depends  strongly  on  pressure.  It  increases  rapidly  as  the  fluid  state  approaches  the 
critical  point,  and  theoretically  becomes  infinite  exactly  at  the  critical  point.  This  implies  that  much 
more  thermal  energy  is  needed  to  heat  up  the  cold  fluid  jet  in  Case  1  (i.e.,  pm  =  A3  MPa )  than  the 
other  two  cases  when  the  fluid  temperature  transits  across  the  near-critical  regime.  The  anomalous 
variations  of  fluid  volumetric  and  thermal  properties  near  the  critical  point  and  their  dependence  on 
pressure  have  profound  influences  on  fluid  jet  development  at  high  pressures. 
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Fig.4  Density  and  constant-pressure  specific  heat  of  nitrogen  as  the  functions  of  temperature 
and  pressure. 


2.6.1  Instantaneous  Flowfield 

For  a  constant-density  jet,  the  shear  layer  between  the  jet  and  the  ambient  fluid  is 
susceptible  to  the  Kelvin-Helmholz  instability  and  experiences  vortex  rolling,  paring,  and  breakup. 
A  cryogenic  supercritical  jet  undergoes  qualitatively  the  same  process,  but  with  additional 
mechanisms  arising  from  volume  dilation  and  baroclinic  torque. 

Figure  5  shows  snapshots  of  the  density,  density-gradient,  temperature,  and  vorticity- 
magnitude  fields  at  three  different  ambient  pressures.  The  small  vortical  structures  in  the  core 
region  result  from  the  imposed  turbulent  motions  at  the  injector  exit.  For  Case  1,  in  which  the 
ambient  pressure  of  4.3  MPa  is  closer  to  the  critical  value,  the  jet  surface  is  straight  near  the  injector 
with  only  tiny  instability  waves  in  the  downstream  region.  As  the  ambient  pressure  increases,  the 
velocity  fluctuation  in  the  radial  direction  becomes  more  vigorous.  Large-scale  instability  waves 
develop  in  the  near-injector  region,  which  then  grow  up  and  roll  into  a  succession  of  ring  vortices  as 
the  injected  fluid  moves  downstream.  The  resultant  vortical  flow  motions  facilitate  the  entrainment 
of  the  ambient  gaseous  nitrogen  into  the  cold  jet  fluid.  The  initial  density-stratification  layer  is  only 
slightly  stretched  in  Case  1,  but  severely  twisted  in  higher-pressure  cases.  As  a  general  trend,  the 
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higher  the  ambient  pressure,  the  stronger  vortical  motions  and  radial  velocity  fluctuations  near  the 
jet  surface. 

It  has  been  established  [28]  that  the  evolution  and  interaction  of  large  coherent  structures 
strongly  influence  the  mixing  and  entrainment  of  a  shear  layer.  Figure  6  shows  the  temporal 
evolution  of  the  jet  surface  structures  of  Case  3  ( />„  =  9.3  MPa ).  Both  the  temperature  and  density 
fields  clearly  demonstrate  the  entrainment  of  lighter  and  warmer  ambient  gaseous  nitrogen  into  the 
jet  flow  through  vortical  motions,  along  with  a  series  of  thread-like  entities  emerging  from  the  jet 
surface.  The  same  phenomena  were  observed  in  the  experiments  by  Chehroudi  and  co-works  under 
the  same  flow  condition  [10]. 


t  1.500  ms  t  -  1.525  ms  t  -  1.550  ms  t  -  1.575  ms 


3.0  6.0  3.0  6.0  3.0  6.0  3.0  6.0 

X-;  Ojnj  -v  X/Djrj 


Fig.6  Time  evolution  of  jet  surface  structures  (pTO  =  9.3  MPa,  TTO  =  300  K,  Ujnj  =15  m/s,  Tinj  = 
120  K,  Dinj  =  254  pm). 


2.6.2  Effect  of  Density  Stratification 

To  explore  the  formation  and  influence  of  the  density-stratification  layer  near  the  jet 
boundary,  conditional-averaged  temperatures  over  the  regions  where  the  density-gradient 
magnitudes  exceed  pre-specified  cutoff  values  are  determined.  The  result  is  listed  in  Table  2,  where 
lv/>L  s is  the  maximum  density  gradient  in  the  entire  field,  around  1.15xl07  kg/m4  for  Case  3. 

The  conditional  fluid  temperature  decreases  with  increasing  cutoff  density  gradient,  and  approaches 
the  inflection  point  on  the  isobaric  p-T  curve  shown  in  Fig.  4.  Table  3  summarizes  the  inflection 
temperatures  for  the  three  different  pressures  considered  herein.  It  is  noteworthy  that  fluid 
properties  usually  undergo  rapid  variations  across  the  temperature  inflection  point  for  an  isobaric 
process,  and  the  specific  heat  reaches  its  maximum  at  this  point.  For  example,  the  density  of 
nitrogen  decreases  more  than  three  times  in  Case  1  as  the  temperature  increases  from  125  to  135  K. 
Thus,  the  formation  of  the  steep  density-gradient  region  is  closely  related  to  property  variations, 
which  to  a  large  extent  are  dictated  by  real-fluid  thermodynamics.  Turbulent  diffusion  and  mixing 
tend  to  introduce  warm  ambient  gases  into  the  cold  jet,  and  subsequently  smooth  the  density- 
stratification  effect.  The  drastic  volume  dilation  during  the  mixing  process  when  the  temperature 
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Table  2  Conditional  averaged  temperatures  (K)  in  regions  with  |Vp|  >  cutoff  value. 


Cutoff 

Case  1 

Case  2 

Case  3 

OHM™ 

131.7 

151.0 

158.0 

0  2 

131.1 

148.0 

154.0 

0-3jVpL 

130.6 

145.5 

151.3 

Table  3  Temperatures  at  inflection  points  on  isobaric  p-T  curves  at  different  pressures 

Case  1  Case  2  Case  3 
Poo  (MPa)  4.3  6.9  9.3 

T  (K)  130.3  142.1  148.6 


transits  across  the  inflection  point,  however,  prevents  the  entrainment  of  the  surrounding  fluid,  and 
thereby  facilitates  the  formation  of  steep  density-gradient  regimes. 

The  effect  of  density  stratification  on  the  evolution  of  a  planar  mixing  layer  was  studied  by 
Atsavapranee  and  Gharib  [29].  Because  the  higher  density  stratification  increasingly  inhibits 
instability-wave  growth  and  vortex  paring,  the  flow  topography  is  considerably  simplified,  as 
evidenced  in  Fig.  5.  In  addition,  the  total  yield  of  mixed  fluid  is  reduced  with  high  density 
stratification  due  to  the  combined  effect  of  weakened  fluid  entrainment  into  the  Kelvin-Helmholtz 
vortices,  decreased  frequency  of  vortex  paring,  and  arrest  of  turbulence  during  flow  restratification. 

Figure  7  shows  the  power  spectral  densities  (PSD)  of  the  axial  and  radial  velocity 
fluctuations  at  various  radial  locations  of  r/Dinj  =  0.6, 0.7,  and  0.8  in  the  mixing  layer  for 

p„  =9.3  MPa .  The  axial  position  of  x/  Dinj  =16  is  near  the  end  of  the  potential  core,  where  a 

sharp  density  gradient  exists  at  r  =  0.lDinj  (see  Fig.  5).  In  the  low-frequency  regime,  in  which 

large-scale  structures  prevail,  the  axial  velocity  fluctuation  increases  as  the  density  stratification 
layer  is  approached.  The  trend  for  the  radial  velocity  fluctuation,  however,  is  opposite.  The  velocity 
fluctuations  in  the  high-frequency  range  remain  basically  insensitive  to  the  radial  position.  Density 
stratification  exerts  a  strong  influence  on  large-scale  flow  motions.  It  acts  like  a  solid  wall  in  the 
flow  that  amplifies  the  axial  turbulent  fluctuation  but  damps  the  radial  one.  The  same  PSD  results 
are  also  observed  for  Cases  1  and  2,  which  are  not  presented  here. 

A  similar  phenomenon  was  reported  by  Hannoun  et  al.  [30]  in  their  experiments  on  grid- 
induced  shear-free  turbulence  near  a  sharp  density  interface.  As  a  consequence  of  the  strong 
anisotropy  of  the  turbulence  near  the  density  interface,  large  eddies  of  integral  length  scales  become 
flattened,  and  the  vertical  component  of  the  turbulent  kinetic  energy  is  transferred  to  its  horizontal 
quantity.  Such  an  energy  redistribution  among  its  spatial  components  considerably  modifies  the 
amount  of  energy  available  for  fluid  mixing  at  the  density  interface.  Thus,  the  existence  of  strong 
density  stratification  suppresses  radial  velocity  fluctuations  in  the  flowfield  and  inhibits  the 
development  of  instability  waves.  In  the  present  study,  the  initial  density  ratio  and  the  strength  of 
density  stratification  decrease  as  the  ambient  pressure  increases,  so  do  their  damping  effects  on  the 
shear  layers.  The  jet  surface  is  nearly  straight  in  Case  1,  with  only  tiny  instability  waves 
developing  in  the  downstream  region.  The  shear  layers,  however,  are  highly  twisted  in  both  Cases 
2  and  3,  as  shown  in  Fig.  5. 
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Fig.7  Power  spectral  densities  of  velocity  fluctuations  at  different  radial  locations  with  x/Di„j  = 
16  (p„o  =  9.3  MPa,  Too  =  300  K,  Ujnj  =15  m/s,  Tinj  =  120  K,  Dinj  =  254  Jim). 


The  frequency  spectra  of  velocity  fluctuations  shown  in  Fig.  7  do  not  indicate  the  standard 
5/3  law  for  the  turbulent  kinetic  energy  spectrum  based  on  the  Kolmogorov  theory,  which 
characterizes  the  inertial  subrange.  Apte  and  Yang  [31]  mentioned  in  their  two-dimensional 
simulations  of  unsteady  flow  evolution  in  a  porous-walled  chamber  that  the  exponent  of  the  wave 
number  in  the  inertial  subrange  of  the  turbulent  energy  spectrum  varies  between  /"3  and  Z-4 ,  with 
/  being  the  frequency.  Gilbert  [32]  proposed  that  the  kinetic  energy  spectrum  can  be  obtained 
from  spiral  vortex  distributions  within  coherent  vortices  and  should  follow  the  /”I1/3  law.  In  this 
work,  the  frequency  dependence  in  the  inertial  subrange  lies  between  /~5/3  and  /“3. 

The  present  two-dimensional  large-eddy  simulation  inherently  neglects  the  vortex  stretching 
mechanism,  which  is  responsible  for  the  transfer  of  energy  from  large  to  small  scales  through 
energy  cascade  and  the  continuous  generation  of  small-scale  vortical  structures.  Since  these 
structures  redistribute  and  dissipate  energy  at  the  smallest  scales,  the  lack  of  vortex  stretching  leads 
to  lower  energy-dissipation  and  turbulence-production  rates  during  the  flow  evolution.  This  issue 
will  be  addressed  in  subsequent  three-dimensional  calculations. 
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2.6.3  Shear-Layer  Instability 

To  study  the  vortical  dynamics  and  flow  instability  in  the  mixing  layer,  the  power  spectral 
densities  of  the  radial  velocity  oscillations  at  two  different  axial  locations  are  presented  in  Fig.  8. 
The  radial  position  is  fixed  at  r/Djnj  =  0.5.  A  dominant  frequency  around  35  kHz,  corresponding  to 
the  most  amplified  frequency  of  the  shear  layer  instability,  is  observed  at  an  upstream  location  of 
x/Dinj  =  8  for  both  Cases  2  and  3.  This  frequency  is  weakly  dependent  on  the  ambient  pressure  (or 
the  density  ratio).  When  the  fluid  is  convected  downstream  to  x/Dinj  =18,  the  dominant  frequency 
decreases  to  15.9  kHz  for  Case  2  and  18.6  kHz  for  Case  3,  nearly  one  half  of  the  value  of  the  most 
amplified  shear  instability  mode.  The  vortex  pairing  process  is  clearly  demonstrated.  The  situation 
with  Case  1,  however,  is  considerably  different.  Owing  to  the  low  compressibility  factor  (i.e.  0.21) 
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Fig.8  Effect  of  pressure  on  power  spectral  densities  of  radial  velocity  fluctuations  at  two 
different  axial  locations  with  r/Dinj  =0.5. 
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of  the  injected  fluid  and  the  large  density  ratio  of  12.24,  the  jet  exhibits  a  liquid-like  fluid  behavior 
distinct  from  the  other  high-pressure  cases.  The  dominant  frequencies  become  20  and  10  kHz  at 
x/Dinj  =  8  and  18,  respectively. 

The  flow  instabilities  and  vortex  shedding  in  constant-density  shear  layers  were  reviewed  by 
Schadow  and  Gutmark  [33].  Based  on  their  work,  the  initial  vortex  shedding  frequency,  fx ,  can  be 
scaled  with  the  shear-layer  momentum  thickness  0O  and  a  characteristic  velocity  U  ,  normally 
taken  as  the  average  bulk  velocity  of  the  two  streams.  The  result  yields  a  non-dimensional 
frequency  or  Strouhal  number,  Sti  =  ft0Q  / U  ,  which  ranges  from  0.044  to  0.048  for  a  planar 
turbulent  shear  layer.  As  the  vortices  move  downstream,  they  merge  together  to  oscillate  at  the  sub¬ 
harmonics  of  the  initial  vortex  shedding  frequency,  fx/  N  (N  =  2,  3,  4,  •  •  •) .  Although  the 

above  analysis  was  formulated  for  planar  flows,  it  can  be  applied  with  good  accuracy  to  mixing 
layers  in  axisymmetric  configurations  if  the  thickness  of  the  shear  layer  is  much  smaller  than  the 
radius  of  the  injector.  In  the  present  work,  U  =8 mis  and  the  initial  momentum  thickness  00 
estimated  for  a  fully  developed  turbulent  pipe  flow  is  0.01 1  mm.  If  the  strouhal  number  is  chosen 
to  be  Stt  =  0.048 ,  then  the  most  amplified  frequency  becomes  /,  =  34.9  kHz  and  the 

corresponding  second  harmonic  frequency  is  17.5  kHz.  Those  values  agree  well  with  the 
calculated  vortex-shedding  frequencies  based  on  the  radial  velocity  oscillations  in  Cases  2  and  3. 

To  provide  more  insight  into  observed  flow  phenomena,  a  linear  stability  analysis  is  carried 
out  of  the  effects  of  ambient  pressure  on  the  fluid  jet  evolution.  The  work  extends  the  approaches 
described  in  Refs.  [34]  and  [35]  to  include  real-fluid  thermodynamics.  The  SRK  equation  of  state  is 
implemented  in  the  formulation.  Each  dependent  variable  is  decomposed  into  a  base  and  a 
perturbation  quantity.  The  former  is  adopted  directly  from  the  present  simulation.  The  latter  takes 
the  following  general  form  for  a  planar  jet. 

<j>(x,y,t)  =  <p(y)zxp{i(kx-ax)}  (15) 


where  k  and  co  are  the  wave  number  and  frequency,  respectively.  For  a  spatial  instability 
problem,  k  is  a  complex  variable  and  its  negative  imaginary  part  represents  the  spatial  growth  rate. 
After  substitution  of  the  decomposed  variables  into  the  conservation  laws  and  linearization  of  the 
result,  a  dispersion  equation  characterizing  the  relationship  between  the  wave  number  and 
frequency  can  be  derived  in  terms  of  pressure  fluctuation  as  follow. 


dy  p  dy  u-colkdy  dy 


(16) 


The  problem  now  becomes  solving  Eq.  (16)  for  the  eigenvalues  k  and  co  subject  to 
appropriate  boundary  conditions.  A  complete  discussion  of  the  stability  analysis  is  given  in  Ref. 
[36]. 


Figure  9  shows  the  spatial  growth  rates  of  the  instability  waves  as  a  function  of  the 
normalized  frequency  (i.e.,  the  Strouhal  number)  for  the  three  different  pressures  considered  in  the 
present  study.  The  growth  rate  of  the  shear  wave  decreases  with  increasing  pressure.  The  strong 
density  stratification  in  the  lower  pressure  case  suppresses  the  growth  of  the  wave  and  stabilizes  the 
mixing  layer.  The  theoretically  predicted  frequencies  of  the  most  unstable  oscillations  are  23.2, 
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28.8  and  29.6  KHz  for  Cases  1,  2,  and  3,  respectively.  These  values  are  slightly  underestimated  for 
Cases  2  and  3,  but  overestimated  for  Case  1.  The  maximum  deviation  from  the  numerical 
simulation  is  15%.  The  frequency  of  the  most  unstable  mode  exhibits  a  weak  pressure  dependence 
at  high  pressures.  It,  however,  decreases  significantly  in  the  near-critical  regime  due  to  the 
enhanced  effect  of  density  stratification  and  increased  mixing-layer  momentum  thickness. 


Fig.9  Spatial  growth  rate  as  function  of  Strouhal  number  at  different  pressures. 


2.6.4  Vortical  Dynamics 

As  a  consequence  of  the  large  velocity  difference  between  the  jet  and  ambient  flow,  a  strong 
shear  layer  is  generated  near  the  injector,  which  is  then  tilted  and  develops  to  large  structures  due  to 
vortex  interactions.  The  vortical  dynamics  can  be  best  quantified  using  the  transport  equation  given 
below. 


—  =  (0V)«-(V-«)fli-V(-)xV/)  +  Vx(-V'r)  (17) 

Dt  pp 

where  D/Dt  stands  for  the  substantial  derivative.  The  first  term  on  the  right-hand  side  represents 
vortex  stretching,  which  vanishes  in  the  present  two-dimensional  simulation.  The  second  term 
describes  the  volume-dilatation  effect,  and  the  third  term  denotes  the  baroclinic  torque  produced  by 
the  misalignment  between  the  pressure  and  density  gradients.  The  last  term  arises  from  viscous 
dissipation.  For  cryogenic  fluid  injection  under  supercritical  conditions,  large  volume  expansion 
occurs  when  the  jet  is  heated  by  the  ambient  gas.  Thus,  both  baroclinic  torque  and  volume 
dilatation  may  play  an  important  role  in  determining  vorticity  transport. 

Figure  10  shows  an  instantaneous  azimuthal  vorticity  budget  in  the  radial  direction  at  the 
axial  position  of  x/  Dinj  =5  for  all  the  three  cases.  At  this  position,  the  large  coherent  structures 

are  well  developed  in  Cases  2  and  3.  The  results  are  normalized  by  the  bulk  velocity  (uinj)  and 

momentum  thickness  ( 0Q )  at  the  injector  exit.  The  baroclinic  torque  and  viscous  dissipation  locally 

rival  in  magnitude,  and  have  opposite  contributions  in  equation  (17).  These  two  terms  attain  their 
maxima  in  regions  with  large  density  gradients  in  which  vigorous  mixing  between  the  jet  and 
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Fig.  10  Vorticity  budgets  at  different  ambient  pressures  (Tco=300  K,  Ujnj=  15  m/s,  Tmj=  120  K,  Djnj 
=  254  pm,  t  =  1.55  ms). 

ambient  flows  occurs  (see  Fig.  5).  In  Case  1,  much  of  the  vorticity  production  takes  place  on  the 
light  fluid  side.  A  similar  phenomenon  was  noted  by  Okong’o  and  Bellan  in  their  study  of  a 
supercritical  binary  mixing  layer.15  As  the  ambient  pressure  increases  (or  the  density  ratio 
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decreases),  the  location  with  intensive  vorticity  production  slightly  shifts  toward  the  dense  fluid,  but 
still  resides  on  the  ambient  gas  side,  since  the  radial  position  of  the  mixing  region  also  changes. 
The  magnitudes  of  all  the  three  vorticity  production  terms  increase  with  increasing  ambient 
pressure  due  to  strengthened  vortical  motions. 

2.6.5  Mean  Flow  Properties 

The  mean  flow  properties  are  obtained  by  taking  long-time  average  of  the  instantaneous 
quantities  over  10  ms  (about  8  flow-through  times)  after  the  calculated  flowfield  has  reached  its 
stationary  state.  Figure  11  shows  the  radial  distributions  of  the  normalized  mean  density, 
p+  =  (p  —  p„)  /(pc  at  six  different  axial  locations.  The  subscript  c  refers  to  the  quantity  at 

the  centerline.  The  radial  coordinate  is  normalized  by  the  full  width  of  the  radial  profile  measured 
where  the  flow  property  of  concern  (i.e.,  density  in  the  present  figure)  is  one  half  of  its  maximum 
value  (FWHM),  ru2 .  Similar  to  incompressible  fluid  jets,  there  exist  three  distinct  regions  in  a 

cryogenic  fluid  jet  under  supercritical  conditions:  a  potential  core  in  the  upstream  where  a  flat-hat 
distribution  around  the  centerline  occurs,  a  transition  region,  and  a  fully  developed  self-similar 
region.  The  density  profiles  at  x/Dinj  >25  merge  to  a  single  distribution  for  Cases  2  and  3, 

manifesting  the  existence  of  self-similarity  in  the  downstream  region.  Such  a  self-similar  profile, 
however,  is  not  observed  until  x/  Dinj  >  30  for  Case  1.  The  high-pressure  condition  facilitates  the 

development  of  the  self-similar  behavior  through  its  effect  on  reducing  the  density  ratio  of  the 
injected  fluid  to  the  ambient  flow.  Figure  12  shows  the  comparison  with  the  density  field  measured 
by  Chehroudi  and  Talley  [12  ]  using  the  Raman  scatting  technique  under  the  same  condition  of 
Case  2  ( =  6.9  MPa ).  The  calculated  and  measured  radial  distributions  follow  exactly  the  same 
trend  at  the  two  positions.  The  maximum  deviation  is  5%. 

Figure  13  shows  the  normalized  axial  velocity  profiles  in  the  radial  direction.  Self¬ 
similarity  is  achieved  at  the  same  location  as  that  for  the  density  distribution.  It  occurs  at  a 
relatively  upstream  position  as  the  ambient  pressure  increases.  Compared  with  the  density  field,  the 
velocity  profile  exhibits  a  narrower  distribution,  with  the  jet  boundary  (defined  as  the  radial 
position  at  which  u  =  0.0  lw(.)  situated  at  rlrU2  =  2.5 ,  a  value  consistent  with  the  data  reported  in 

Ref.  [9].  In  addition,  no  flat-hat  profile  exists  even  in  the  upstream  region  due  to  the  use  of  the  one- 
seventh  power  distribution  for  a  fully  developed  turbulent  pipe  flow  at  the  injector  exit. 

Figure  14  presents  the  axial  distributions  of  the  normalized  temperature,  density  and 
compressibility  factor  along  the  centerline  for  three  different  pressures.  Evidently,  the  temperature 
decreases  relatively  slowly  with  decreasing  pressure  in  the  axial  direction.  It  is  well  established  that 
when  a  fluid  reaches  its  thermodynamic  critical  state,  the  constant-pressure  specific  heat  becomes 
infinite  and  the  thermal  diffusivity  decreases  to  zero,  a  phenomenon  known  as  the  critical 
divergence.  Because  the  ambient  pressure  for  Case  1  ( px  =  4.3  MPa )  is  closer  to  the  critical  value, 
the  specific  heat  increases  drastically  when  the  temperature  transits  across  the  inflection  point  on 
the  isobaric  p-T  curve  (see  Fig.  4).  This  effect,  combined  with  the  lower  thermal  conductivity, 
causes  a  slower  increase  in  temperature  along  the  centerline  in  Case  1  as  the  injected  fluid  moves 
downstream.  The  compressibility  factor,  shown  in  Fig.  2,  indicates  a  monotonic  decrease  with 
decreasing  pressure  in  the  low-temperature  regime,  and  a  rapid  increase  across  the  inflection  point 
on  the  isobaric  p-T  curve.  The  compressibility  factor,  thus,  increases  much  more  rapidly  at  a  lower 
pressure,  which  amounts  to  a  faster  decrease  in  density.  Although  the  centerline  temperature  varies 
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Fig.  1 1  Effect  of  pressure  on  normalized  density  distribution  in  radial  direction. 


laggardly  in  all  the  three  cases,  the  rapid  variation  of  the  compressibility  factor  causes  the  density  to 
undergo  a  fast  decrease  downstream  of  the  potential  core,  especially  for  the  case  of  p„  =  4.3  MPa . 

The  rapid  variations  of  thermophysical  properties  exert  a  significant  influcence  on  the  jet 
flow  evolution.  Figure  15  shows  the  radial  distributions  of  the  specific  heat,  thermal  diffusivity, 
and  kinematic  viscosity  at  x!  Dinj  =10,  a  location  slightly  upstream  of  the  end  of  the  potential  core. 
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Fig.  12  Comparison  of  calculated  and  measured  density  profiles  at  different  axial  locations. 

A  spike  exists  in  the  specific-heat  profile  at  rl  rll2  =  0.6  in  Case  1,  mainly  due  to  the  temperature 

transition  across  the  inflection  point.  The  same  phenomenon  is  observed  for  Cases  2  and  3,  but 
with  much  lower  amplitudes.  The  thermal  diffusivity  also  exercises  a  large  excursion  of  variation 
across  the  temperature  inflection  point.  It  increases  considerably  from  the  liquid  core  to  the 
ambient  flow  by  a  factor  of  14,  7,  and  5  in  Cases  1,  2,  and  3,  respectively.  A  similar  situation 
occurs  with  the  kinematic  viscosity.  The  intensive  change  of  viscosity  in  the  shear  layer  gives  rise 
to  another  vorticity  production  mechanism,  as  shown  in  equation  (17).  These  observations,  again, 
highlight  the  importance  of  thermophysical  properties  in  dictating  the  behavior  of  a  supercritical 
fluid  jet  [37], 

2.7  Conclusions 

A  comprehensive  numerical  study  has  been  conducted  to  investigate  cryogenic  fluid 
injection  and  mixing  under  supercritical  conditions.  The  model  accommodates  full  conservation 
laws  and  real-fluid  thermodynamics  and  transport  phenomena  over  the  entire  range  of  fluid  states  of 
concern.  Turbulent  closure  is  achieved  using  a  large-eddy-simulation  technique.  The  present 
analysis  allows  a  detailed  investigation  into  the  temporal  and  spatial  evolution  of  a  cryogenic  jet. 
The  near-field  behavior'  is  well  captured. 

The  major  results  obtained  are  summarized  below. 

1.  As  a  result  of  intensive  property  variations  between  the  fluid  jet  and  surroundings,  a  series 
of  large  density-gradient  regions  are  formed  around  the  jet  surface.  These  regions  act  like  a 
solid  wall  that  amplifies  the  axial  flow  oscillations  but  damps  the  radial  ones.  The  interfacial 
instability  in  the  shear  layer  is  effectively  suppressed,  especially  for  cases  with  large  density 
ratios.  As  the  ambient  pressure  increases,  the  strength  of  density  stratification  decreases,  so 
does  its  damping  effect.  Thus,  the  jet  expands  rapidly  with  increasing  pressure. 

2.  Various  mechanisms  dictating  vorticity  transport  are  analyzed.  The  baroclinic  torque 
arising  from  the  density  stratification  between  the  injected  and  ambient  flows  and  viscous 
dissipation  play  an  important  role  in  determining  the  flow  evolution. 


3.  The  jet  dynamics  are  largely  dictated  by  the  local  thermodynamic  state  of  the  fluid.  When 
the  temperature  transits  across  the  inflection  point  in  an  isobaric  process,  the  rapid  property 
variations  may  qualitatively  change  the  jet  behavior  compared  with  its  counterpart  at  low 
pressures.  In  addition,  an  increase  in  the  ambient  pressure  may  result  in  an  earlier  transition 
of  the  jet  into  the  self-similar  region. 

4.  The  spatial  growth  rate  of  the  surface  instability  wave  increases  as  the  ambient  pressure 
increases.  The  frequency  of  the  most  unstable  mode  exhibits  a  weak  pressure  dependence  at 
high  pressures.  It,  however,  decreases  significantly  in  the  near-critical  regime  due  to  the 
enhanced  effect  of  density  stratification  and  increased  mixing-layer  momentum  thickness. 
The  result  agrees  well  with  the  linear  stability  analysis. 
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Fig.  15  Radial  distributions  of  mean  thermophysical  properties  at  x/Djnj  =  10. 
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Task  3 

Shear  Co-axial  Injection  of  Cryogenic  Fluids  with  Acoustic  Excitations 

This  task  simulated  the  experiments  conducted  by  Chehroudi,  Talley,  and  colleagues  at  the 
Air  Force  Research  Laboratory  (AFRL).  The  purpose  is  to  investigate  cryogenic  fluid  injection  and 
mixing  of  a  shear  coaxial  injector  under  conditions  with  and  without  externally  impressed  acoustic 
forcing.  A  cryogenic  nitrogen  jet  is  injected  through  a  sharp-edged  stainless-steel  tube  having  a 
length  of  50  mm  and  an  inner  diameter  0.508  mm.  The  inner  tube  thickness  is  0.415  mm.  The 
outer  diameter  of  the  annular  stream  is  2.42  mm.  To  facilitate  the  application  of  acoustic  boundary 
condition  correctly,  a  cuboid  computational  domain  is  generated.  The  entire  grid  system  (except  for 
the  region  inside  the  injector)  has  141x81x301  (3.44  million)  cells  along  the  axial,  radial  and 
azimuthal  directions,  respectively,  of  which  41  axial  cells  are  used  to  cover  the  LN2  inlet  section 
and  31  axial  cells  for  the  GN2  inlet  section.  Square  grids  are  implemented  in  the  center  region  in 
order  to  avoid  the  numerical  singularity  at  the  centerline  of  the  computation  domain  and  to  increase 
the  computational  efficiency.  The  grid  resolution  is  chosen  based  on  the  inlet  Reynolds  number, 
such  that  the  largest  grid  size  falls  in  the  inertial  sub-range  of  the  turbulent  energy  spectrum.  The 
analysis  is  conducted  on  an  in-house  distributed-memory  parallel  computer.  The  computational 
domain  is  divided  into  72  blocks,  and  a  total  number  of  57  processors  are  used. 

Results  from  this  task  have  led  to  the  following  technical  publication. 

“Dynamics  of  Shear-Coaxial  Cryogenic  Nitrogen  Jets  with  Acoustic  Excitation  under 

Supercritical  Conditions,”  by  T.  Liu,  N.  Zong,  and  V.  Yang,  AIAA  Paper  2006-0759, 

presented  at  the  44th  AIAA  Aerospace  Sciences  Meeting  and  Exhibit,  January  2006. 

3.1  Introduction 

Combustion  instability  has  long  been  one  of  the  most  complicated  phenomena  in  the 
development  of  liquid  rocket  engines.  The  prevalence  of  instabilities  is  primarily  attributed  to  two 
fundamental  phenomena^  1]  Combustion  chambers  are  almost  entirely  closed  and  the  internal 
processes  tending  to  attenuate  unsteady  motions  are  weak  and  the  energy  required  to  drive  unsteady 
motions  represents  exceedingly  small  fraction  of  the  heat  released  by  combustion.  Characterized  by 
high  amplitudes  and  high  frequencies,  acoustic  instability  is  the  most  destructive.  It  will  cause  local 
burnout  of  combustion  chamber  walls  and  injector  plates  and  may  significantly  shorten  the  lifetimes 
of  a  combustor[2].  In  order  to  thoroughly  understand  the  nature  of  interactions  between  acoustic 
waves  and  liquid  propellant  combustion  in  rocket  engines,  it  is  important  to  conduct  theoretical  and 
experimental  studies  on  simplified  configurations. 

Injectors  play  an  important  role  in  defining  the  flow  structure  and  flame  dynamics.  In 
modem  liquid  rocket  engines,  shear  coaxial  jet  injectors  have  been  commonly  used.  Coaxial 
atomization  can  be  divided  into  five  distinctive  physical  regimes[3]:  1)  turbulent  liquid  core;  2) 
liquid/gas  interface  where  primary  atomization  occurs;  3)  dense  spray  where  second  atomization 
and  droplet  breakup  occur;  4)  dilute  spray  where  droplets  are  dispersed;  5)  gas  flow  region.  The 
dominant  flow  parameters  involved  include:  Reynolds  number  (Re),  momentum  flux  ratio  (M),  area 
ratio,  and  Weber  number  (We).  The  primary  advantage  of  a  shear  coaxial  jet  injector  lies  in  the  fact 
that,  with  a  large  outer(fuel)-to-inner(oxidizer)  momentum  flux  ratio,  it  is  efficient  to  fragment  and 
mix  the  reactants  to  a  satisfactory  level  of  uniformity  for  evaporation  and  combustion  within  a 
preferred  distance. 
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Studies  of  shear  coaxial  injection  have  been  extensive.  Mayer  and  Kriille[3]  found  in  their 
experiment  on  water  and  air  that  similar  to  increasing  gas  velocity,  increasing  gas  density  results  in 
a  faster  and  finer  atomization.  The  phenomena  may  be  attributed  to  the  increased  surface  wave 
growth  and  instabilities  of  liquid  jets  and  enhanced  liquid  atomization.  Rehab,  Villermaux,  and 
Hopfinger[4]  studied  the  effects  of  large  velocity  ratio  on  coaxial  water/water  jets.  When  the 
velocity  ratio  of  the  annular  to  the  inner  jet  is  less  than  a  critical  value,  the  inner  potential  core  is 
found  to  be  inversely  dependent  on  the  velocity  ratio.  Otherwise,  the  inner  potential  core  will  be 
truncated  by  the  unsteady  recirculating  bubble.  The  same  authors[5]  later  examined  the  geometric 
effects.  A  larger  annular  gap  was  found  to  result  in  a  longer  inner  potential  core  and  a  larger  value 
of  the  critical  velocity  ratio.  The  location  of  the  recirculation  zone  appeared  at  a  further  location. 
An  increase  of  the  retraction  of  the  inner  jet  exit  increases  inner  potential  core,  and  the  recirculation 
zone  occurs  closer  to  the  outer-jet  exit.  Lasheras,  Villermaux,  and  Hopfinger[6]  conducted 
experiments  on  water  jets  surrounded  with  co-flowing  high-speed  annular  air.  Weber  number  was 
found  to  have  a  weak  effect  on  the  break-up  when  We  >  200.  The  momentum-flux  ratio  exerted  a 
dominant  effect  on  the  liquid-intact  length,  especially  with  a  larger  initial  Weber  number. 

So  far,  limited  effort  was  devoted  to  the  effect  of  acoustic  excitation.  Chehroudi,  Davis  and 
co-workers[7,8,9]  studied  the  interactions  between  acoustic  waves  and  a  cryogenic  jet  at  sub-  and 
supercritical  pressures.  The  effects  of  acoustic  waves  on  the  coaxial  injection  with  liquid  and 
gaseous  nitrogen  were  examined  over  a  wide  range  of  chamber  pressures.  It  was  found  that  the 
acoustic  effects  were  substantial  under  subcritical  conditions,  but  became  unnoticeable  as  pressure 
increased  above  the  critical  point.  The  shear  coaxial  jets  at  subcritical  conditions  were  more  easily 
to  be  influenced  with  a  higher  velocity  ratio. 

The  present  work  attempts  to  simulate  the  experiment  work  described  in  Ref.7-9.  It  deals 
with  numerical  simulation  of  coaxial  fluid  injection  and  mixing  dynamics  under  super-critical 
conditions  with/without  acoustic  excitation.  Results  will  enhance  basic  understanding  of 
interactions  between  acoustic  waves  and  coaxial  jets  in  cryogenic-propellant  rocket  engines.  This 
portion  of  work  is  organized  as  follows.  Section  II  and  III  describe  the  theoretical  formulation  and 
numerical  method,  respectively.  In  Section  IV,  the  computational  domain  and  grid  system  is 
addressed,  and  primary  results  are  presented  and  discussed.  Finally,  a  brief  summary  is  given  in 
Section  V. 

3,2  Theoretical  Formulation 

Several  computational  and  modeling  challenges  must  be  overcome  in  order  to  perform 
numerical  simulations  on  high-pressure  cryogenic  fluid  dynamics.  First,  thermodynamic  non¬ 
idealities  and  transport  anomalies  occur  as  the  fluid  transits  through  the  transcritical  regime.  Thus 
treating  these  phenomena  in  a  manner  consistent  with  the  intrinsic  characteristics  of  a  numerical 
algorithm  presents  a  main  obstacle.  Second,  the  rapid  variation  of  fluid  state  and  wide  disparities  in 
the  characteristic  time  and  length  scales  impose  the  well-known  stiffness  problem.  In  order  to 
circumvent  all  these  difficulties,  a  comprehensive  theoretical  framework  capable  of  treating 
turbulence,  transcritical  property  variations,  and  general  fluid  thermodynamics  is  developed. 

3.2.1  Governing  Equations 

The  present  analysis  is  based  on  a  large-eddy-simulation  (LES)  technique,  in  which  large- 
scale  turbulent  structures  are  directly  computed  and  small  dissipative  structures  are  modeled. 
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Mathematically,  the  LES  methodology  begins  with  filtering  of  small-scale  effects  from  large-scale 
motions  in  the  full  conservation  equations. 


/  =  /  +  /' 


(1) 


The  Favre-filtered  averaging  is  employed  here  to  simplify  the  governing  equations  and  to  account 
for  the  variable  density  effects. 


(2) 


The  Favre-filtered  conservation  equations  of  mass,  momentum,  and  energy  can  be  expressed 
in  the  following  conservation  form: 
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where  p,  m„  p,  eh  qj  and  %  represent  the  density,  velocity  components,  pressure,  specific  total 
energy,  heat  flux,  and  viscous  stress  tensor,  respectively.  The  subgrid-scale  terms  in  Eqs.  (4)-(5), 
i.e.,  the  subgrid  stress,  subgrid  energy  fluxes,  and  subgrid  viscous  work  are  given  below. 
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The  sgs  terms  are  closed  by  implementing  an  improved  Smagorinsky  model  proposed  by 
Erlebacher  et  al.[10] 
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where 
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n=C,A2|5| 
k"  =  C,  A25,y5.. 


(11) 


(12) 

The  dimensionless  quantities  C«  and  C/  represent  the  compressible  Smagorinsky  constants,  A  is  the 
filter  width,  St]  =  (dui/dxj  +  duj/dxi}/ 2 ,  and  Prf  the  turbulent  Prandtl  number  for  which  a  standard 

value  of  0.7  is  used.  The  sgs  viscous  work  term  crT ,  is  neglected  in  the  present  work  due  to  its 
small  contribution  to  the  total  energy  equation. 

3.2.2  Equation  of  State  and  Thermodynamic  Properties 

A  modified  Soave-Redlich-Kwong  (SRK)  equation  of  state  is  selected  in  the  present  study 
because  of  its  reasonable  accuracy  over  the  high-pressure  and  low-temperature  regime  in  the 
present  study  and  ease  of  implementation!  1 1], 

RT  aa 

P~(* 7Zbj~V(v7bj  (13) 

where  V  and  R  are  the  molar  volume  and  universal  gas  constant,  respectively.  The  energy  and 
volume  constants,  a  and  b,  are  determined  respectively  from  the  following  universal  relations  in 
terms  of  the  critical  temperature  and  pressure. 

a  =  0.42747 R%2  /  Pc \b  =  0.0&664RTJ  pc  (14) 

The  parameter  or  is  a  function  of  the  reduced  temperature  and  acentric  factor. 

a  =  (l  +  (0.48508  +  1.55171ft>-0.15613(y2)(l  -7)05))2 ;  Tr  =  T/Tc  (15) 

In  order  to  avoid  computationally  intensive  iterations  in  calculating  thermodynamic 
properties,  an  accurate  curve  fit  of  T  (e,p)  is  employed  to  update  the  temperature  based  on  the 
specific  internal  energy  and  density  over  the  entire  fluid  state  of  concern.  Thermodynamic 
properties  including  enthalpy,  internal  energy  and  heat  capacity  are  directly  calculated  by  means  of 
fundamental  thermodynamic  theories. 

3.2.3  Transport  Properties  Evaluation 

Transport  properties  including  viscosity  and  thermal  conductivity  are  estimated  by  means  of 
the  corresponding  state  principles  along  with  the  use  of  the  32-term  Benedict-Webb-Rubin  (BWR) 
equation  of  state,  which  are  also  implemented  using  curve  fits  of  ji  (p,  T)  and  A.  (p,  T). 

3.3  Numerical  Method  and  Boundary  Conditions 

The  theoretical  formulation  outlined  above  is  solved  by  means  of  a  density-based,  finite 
volume  methodology.  The  spatial  discretization  employs  a  fourth-order,  central-difference  scheme 
in  generalized  coordinates.  A  fourth-order  scalar  dissipation  with  a  total-variation-diminishing 
switch  developed  by  Swanson  and  Turkel  is  implemented  to  ensure  computational  stability  and  to 
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prevent  numerical  oscillations  in  regions  with  steep  gradients.  Temporal  discretization  is  obtained 
using  a  four-step  fourth-order  Runge-Kutta  scheme.  A  multi-block  domain  decomposition 
technique,  along  with  static  load  balance,  is  employed  to  facilitate  the  implementation  of  parallel 
computation  with  message  passing  interfaces  at  the  domain  boundaries. 

The  boundary  conditions  used  in  the  present  work  are  listed  below.  At  the  inlet  boundary, 
the  three  velocity  components  and  temperature  are  specified.  The  pressure  is  obtained  from  a 
simplified  one-dimensional  axial-momentum  equation.  At  the  outflow  exit,  the  ambient  pressure  is 
specified  and  the  three  velocity  components  and  temperature  are  extrapolated  from  the  interior.  The 
non-slip,  adiabatic  wall  conditions  are  applied  to  all  solid  walls.  In  the  case  with  acoustic 
excitations,  periodic  pressure  oscillations  are  imposed  at  the  side  wall  in  order  to  generate 
transverse  standing  acoustic  waves  in  the  chamber. 

3.4  Results  and  Discussions 

3.4.1  Computational  Case  Description 

The  analysis  simulates  the  experiments  by  Chehroudi,  Davis  and  coworkers  [7-9].  Liquid 
nitrogen  is  delivered  through  the  inner  tube  and  gaseous  nitrogen  flows  through  the  concentric  tube 
at  a  higher  velocity.  The  temperatures  of  inner  and  outer  jets  are  specified  as  132  K  and  191  K, 
respectively.  The  diameters  of  the  inner  and  outer  tubes  are  0.508  mm  and  1.59  mm,  respectively 
and  the  annulus  height  is  0.415  mm  as  shown  in  Fig.  1. 


Figure  1.  Schematic  of  the  coaxial  injector  employed  for  the  high-pressure  LN2/GN2  mixing 
simulation.  (R1  =  0.254mm,  R2  =  0.795mm,  R3  =  1.21mm) 

In  order  to  efficiently  implement  acoustic  boundary  conditions,  a  cuboid  computational 
domain  shown  in  Fig.  2  is  employed.  For  cases  without  acoustic  excitation,  the  main  body  grid 
system  besides  the  inlet  and  buffer  parts  has  301x141x81  (3.44  million)  points  along  the  axial, 
radial  and  azimuthal  direction,  respectively,  of  which  3 1  radial  points  are  used  to  cover  the  LN2  and 
GN2  inlet  section,  respectively.  Square  grids  are  used  in  the  center  part  to  avoid  the  numerical 
singularity  in  the  center  region.  The  grid  resolution  is  chosen  based  on  the  inlet  Reynolds  number 
such  that  the  largest  grid  size  falls  in  the  inertial  sub-range  of  the  turbulent  energy  spectrum.  All 
the  calculations  are  performed  on  an  in-house  distributed-memory  parallel  computer.  The 
computation  domain  is  divided  into  72  blocks,  and  a  total  number  of  57  processes  are  used.  For 
cases  involving  acoustic  excitations,  the  computational  domain  is  extended  in  the  y  direction  in 
order  to  generate  standing  acoustic  waves  with  desired  frequencies. 
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Acoustic  Boundary 
Conditions 


Units:  mm 


Figure  2.  Schematic  of  grid  system,  total  grid  points,  141x81x301  =  3.44  million.  (The  presented 
grid  has  fewer  points  than  these  used  in  the  calculations,  but  the  distributions  of  grid 
points  are  similar) 

Four  different  flow  conditions  without  acoustic  excitation  were  first  treated  to  investigate 
the  effects  of  ambient  pressure,  velocity  ratio  and  momentum  flux  ratio  on  the  injector  flow 
evolution.  Table  1  summarizes  the  flow  parameters  and  simulation  conditions.  The  temperature 
selected  here  matches  the  supercritical  cases  studied  by  Chehroudi,  et  al.  [7-9].  Case  1  corresponds 
to  one  of  their  experimental  conditions  with  a  chamber  pressure  of  4.94  MPa.  The  velocities  of  the 
inner  and  out  jets  are  increased  proportionally  to  shorten  the  computational  time  while  remaining 
the  same  momentum-flux  and  velocity  ratios  as  those  in  the  experiments.  The  chamber  pressure  of 
case  2  is  increased  to  10  MPa  with  all  the  others  conditions  remaining  fixed  in  order  to  examine  the 
effect  of  chamber  pressure.  Cases  3  and  4  have  the  same  conditions  as  cases  1  and  2,  respectively, 
except  that  different  annular  flow  velocities  are  selected  in  order  to  obtain  different  momentum  flux 
ratios.  In  the  end,  cases  1  and  4  have  momentum-flux  ratios  about  3.45,  same  as  that  of  the  SSME 
prebumer  injector.  Case  2  has  a  momentum  flux  ratio  1.0,  which  is  close  to  that  of  the  SSME  main 
chamber  injector. 

A  simple  2D  computational  domain,  as  shown  in  Fig.  3,  along  with  periodic  boundary 
conditions  in  the  z  direction  was  first  employed  to  test  the  specification  of  acoustic  boundary 
conditions.  The  chamber  pressure  is  set  as  4.94  MPa  and  temperature  as  233  K.  Sinusoidal 
pressure  oscillations  with  a  frequency  3000Hz  are  imposed  at  one  of  the  side  wall.  Figure  4  shows 
the  generated  acoustic  velocity  probed  at  different  locations  along  line  2  at  different  times.  Figure  5 
shows  the  acoustic  velocity  and  pressure  along  the  y  direction.  It  can  be  concluded  from  both 
figures  that  the  acoustic  waves  generated  are  nearly  perfect  standing  waves.  The  deviation  may 
arise  from  the  real-fluid  effects,  which  leads  to  a  more  complex  relationship  between  fluid  density 
and  temperature.  Selected  frequency  of  3000  Hz  is  exactly  not  exactly  natural  frequency  of  the 
chamber  in  the  y  direction  because  of  the  presence  of  real  fluid  density  in  the  inner  jet  region.  The 
largest  velocity  oscillation  occurs  at  the  acoustic  antinodal  line  in  the  inner  jet  regions. 


59 


Table  1  Simulation  conditions  for  analysis  of  shear-coaxial  cryogenic  nitrogen  mixing  process. 


case  1 

case  2 

case  3 

case  4 

P  (MPa) 

4.94 

10 

4.94 

10 

Tchm(  K) 

233 

233 

233 

233 

Tf( K) 

191 

191 

191 

191 

T0  (K) 

132 

132 

132 

132 

#(kg/m3) 

98.8 

217.2 

98.8 

217.2 

p0  (kg/m3) 

404.0 

555.4 

404.0 

555.4 

«/(m/s) 

120 

120 

65 

95 

u0  (m/s) 

32 

32 

32 

32 

«//  u0 

3.75 

3.75 

2.03 

2.97 

(pu)fl  (pu)0 

0.92 

1.47 

0.50 

1.16 

(pu2)f  / 

(pu2)o 

3.44 

5.50 

1.01 

3.45 

11.9 

19.0 

6.4 

15.0 

a/  (m/s) 

279.6 

302.2 

279.6 

302.2 

a„  (m/s) 

230.9 

441.2 

230.9 

441.2 

Mf 

0.43 

0.40 

0.23 

0.31 

M0 

0.14 

0.08 

0.14 

0.08 

Re/ 

3.3E5 

5.8E5 

1.8E5 

4.6E5 

Re„ 

1.3E5 

1.0E5 

1.3E5 

1.0E5 

Re  based  on  /?,  and  R0 

Figure  3.  Schematic  diagram  of  acoustic  generation 
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cycle 


Figure  4.  Normalized  acoustic  velocities  at  different  locations  on  line  2. 
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Figure  5.  Instantaneous  acoustic  velocity  and  pressure  along  probe  line  2. 


3.4.2  Flow  Evolution  without  External  Acoustic  Excitation 

Investigations  were  first  conducted  into  injector  dynamics  under  conditions  without 
externally  imposed  acoustic  oscillations.  Figure  6  shows  the  snapshots  of  the  temperature  and 
temperature-gradient  fields.  The  near-field  injector  flow  dynamics  can  be  characterized  by  the 
evolution  of  three  mixing  layers.  The  inner  mixing  layer  forms  around  the  liquid  nitrogen  jet  and 
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Figure  6.  Effect  of  pressure  and  velocity  ratios  on  temperature,  and  temperature- gradient  fields. 
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the  other  two  emerge  from  the  inner  and  outer  boundaries  of  the  annular  gaseous  nitrogen  stream. 
Two  shear  layers  separated  by  the  port  merge  together  after  large-scale  vortices  originating  from  the 
outer  rim  of  the  tube  grow  up.  Similar  to  experimental  observations,  the  evolution  of  the  outer 
shear  layers  prevails  over  the  inner  one  because  of  the  greater  velocity  difference  between  the  high¬ 
speed  outer  stream  and  the  surroundings.  At  the  same  chamber  pressure,  (i.e.  cases  1  and  3),  the 
inner  cold  fluid,  which  is  colored  by  blue,  fades  away  after  a  short  flow  path  for  the  case  with  a 
faster  annular  stream  (case  1).  The  phenomenon  may  be  attributed  to  the  enhanced  turbulent 
mixing  at  a  higher  velocity  ratio.  At  a  higher  velocity  ratio,  the  large-scale  vortices  pinch  the 
centerline  earlier  and  the  potential  cores  of  the  inner  streams  are  reduced.  Figure  7  shows  the 
mean-axial  velocity  and  temperature  fields.  For  all  cases,  the  existence  of  an  inner  and  an  outer 
potential  core  are  clearly  demonstrated.  The  longest  inner  potential  core  is  obtained  in  case3,  owing 
to  the  strong  density  stratification  and  the  lower  momentum  flux  ratio  between  the  inner  and  outer 
streams.  Jet  spreading  angle  is  an  important  parameter  to  measure  the  mixing  and  entrainment 
between  the  injected  flow  and  the  ambience.  At  a  constant  chamber  pressure,  the  outer  jet 
spreading  angle  increases  with  increasing  of  annular  jet  velocity.  The  existence  of 
wake/recirculation  flow  region  immediately  behind  the  post  is  manifested  by  the  negative  velocity 
contour  in  Fig.  7. 


Figure  7.  Effect  of  momentum  flux  ratio  on  flow  evolution,  illustrated  by  contours  of  mean  axial 
velocity  and  temperature  for  different  cases. 
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Figure  8  shows  the  streamlines  based  on  the  mean  axial  and  radial  velocities.  Flow 
recirculation  clearly  occurs  between  the  inner  and  annular  jets.  Rehab,  Villermaux  and  Hopfinger 
[4]  discussed  the  mechanism  responsible  for  the  occurrence  of  the  recirculating  regime  on  the  jet 

axis.  When  the  jet  velocity  ratio  is  large  enough,  the  dynamic  pressure  pu'  /2  is  not  sufficiently 
strong  to  accelerate  the  flow  in  the  positive  direction  and  a  reverse  flow  thus  occurs.  A  model  was 
also  proposed  to  calculate  the  critical  velocity  ratio.  The  velocity  ratios  in  this  work,  however,  are 
much  less  than  their  predictions  because  of  real-fluid  effects  and  geometrical  difference.  In  their 
investigation,  the  thickness  of  the  inner  tube  is  assumed  to  be  zero  and  the  recirculation  bubble 
appears  on  the  jet  axis.  In  the  present  analysis,  the  thickness  of  the  inner  tube  is  finite.  The 
resultant  wake  can  easily  give  rise  to  flow  reversal  behind  inner  tube.  The  same  phenomena  were 
found  by  Davis[12].  He  described  a  recirculation  zone  attached  to  the  lip  of  the  inner  tube,  which 
agrees  well  with  our  simulation.  If  we  further  increase  the  momentum-flux  ratio  to  a  sufficiently 
large  value,  recirculation  zone  may  extend  to  the  jet  axis  line. 


Figure  8.  Mean  axial  velocity  field  and  streamlines  based  on  mean  axial  and  radial 
components. 


velocity 


Figure  9  shows  snapshots  of  the  isosurfaces  of  positive  Q  and  pressure  for  case  3.  Positive 
Q  isosurfaces  isolates  areas  where  the  strength  of  rotation  overcomes  the  strain,  thus  making  those 
surfaces  eligible  as  vortex  envelopes!  13].  It  is  defined  as 
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Large-scale  helical  structures  and  discrete  vortex  rings  form  intermittently  in  the  near  field  from  the 
outer  shear  layer.  These  rings  seem  to  have  the  same  wavelength  until  they  begin  to  disappear  in 
the  downstream  region.  Figure  9(a)  is  plotted  with/without  showing  parts  outside  the  inner  tube. 
The  same  structures  are  devoted  by  the  sky-blue  color,  representing  the  locations  of  recirculation 
bubbles.  Figure  10  shows  the  power  spectral  densities  (PSD)  of  the  pressure  oscillation  at  three 
near-field  locations  in  inner,  annular  jets  and  recirculation  region,  respectively.  The  result  suggests 
the  oscillation  frequency  of  the  recirculation  region  is  determined  by  the  vortex  shedding  frequency. 
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Figure  9.  Large  coherent  structures  visualized  by  iso-surface  of  instantaneous  Q  (a)  and  pressure 
(b)  at  4.935MPa  for  case  3.  (color  contours  are  instantaneous  temperature) 


3.4.3  Flow  Evolution  with  External  Acoustic  Excitation 

To  study  the  effects  of  acoustic  waves,  external  pressure  oscillation  with  an  amplitude 
15kPa  and  a  frequency  3000Hz  was  imposed  at  one  sidewall  after  the  flowfield  becomes  stationary. 
The  power  spectral  densities  of  the  pressure  and  velocity  oscillation  near  the  opposite  sidewall  were 
conducted.  The  resultant  dominant  frequency  of  3000Hz  is  identical  to  the  frequency  of  externally 
imposed  pressure  oscillations,  confirming  the  success  of  acoustic  generation. 

Figure  1 1  shows  the  evolution  of  the  temperature  field  within  one  cycle  of  oscillation  (T  ~ 
0.33  ms)  in  the  slices  perpendicular  to  the  acoustic  velocity  direction.  Even  an  acoustic  excitation 
with  very  low  amplitude  could  substantially  modify  the  injector  flow  dynamics.  Some  sinuous 
wave-like  structures  are  clearly  observed  on  the  jet  surface.  Figure  12  shows  the  evolution  of  the 
temperature  field  in  the  slices  parallel  to  the  acoustic  velocity  direction.  No  obvious  sinuous 
structure  exists;  suggesting  that  oscillatory  motion  of  the  jet  is  primarily  in  the  direction 
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Figure  11. Consecutive  snapshots  of  temperature  field  in  the  slices  perpendicular  to  the  acoustic 
velocity  direction.  Time  increases  with  an  interval  of  57.6|is  between  frames. 
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t-  6.916  ms 


t  =  6.974  ms 


Figure  12.  Consecutive  snapshots  of  temperature  field  in  the  slices  parallel  to  the  acoustic  velocity 
direction.  Time  increases  with  an  interval  of  57.6ps  between  frames. 

3.5  Conclusion 

A  comprehensive  numerical  analysis  has  been  conducted  to  investigate  shear  coaxial 
injection  and  mixing  of  cryogenic  and  gaseous  nitrogen  under  supercritical  pressures.  The  model 
accommodates  full  conservation  laws,  real-fluid  thermodynamics  and  transport  phenomena  over  the 
entire  fluid  state  of  concerned.  The  influences  of  ambient  pressure  and  momentum  flux  ratio  on  the 
injector  flow  dynamics  are  systematically  investigated.  The  effects  of  external  acoustic  forcing  on 
the  flow  evolution  are  also  carefully  examined.  Various  important  flow  features  that  greatly 
influenced  the  injector  dynamics  are  captured.  Major  results  are  summarized  below. 

1.  The  near-field  injector  flow  dynamics  could  be  characterized  by  the  evolution  of  three  shear 
layers  originating  from  the  rims  of  the  two  concentric  tubes.  Due  to  the  greater  velocity 
difference  between  the  high-speed  outer  stream  and  the  surroundings,  the  evolution  of  the  outer 
shear  layers  prevails  over  the  inner  one. 

2.  A  wake/recirculation  region  is  observed  immediately  behind  the  post.  It  oscillates  at  a  frequency 
the  same  as  the  dominant  vortex  shedding  frequency  of  the  shear  layer.  The  interaction  between 
the  wake  region  and  the  annular  jet  is  revealed. 

3.  Owing  to  the  intensive  fluid  property  variations  in  the  flow  field,  a  series  of  large  density- 
gradient  regions  are  formed  around  the  jet  surface.  The  existence  of  those  large  density- 
gradient  regions  significantly  influences  the  injector  flow  dynamics. 

4.  The  dominate  effects  of  the  momentum  flux  ratio/velocity  ratio  on  the  flow  evolution  are 
demonstrated.  As  the  injection  velocity  of  the  annular  stream  increases,  turbulent  mixing  is 
greatly  enhanced.  The  inner  and  outer  potential  cores  are  reduced. 
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5.  The  influence  of  ambient  pressure  on  the  jet  mixing  and  entrainment  is  also  examined.  An 
increase  of  chamber  pressure  gives  rise  to  a  shorter  inner  potential  core  and  a  smaller  jet 
spreading  angle.  The  trend  is  consistent  with  the  experimental  observations. 

6.  The  external  acoustic  excitation  has  substantial  influence  on  the  injector  flow  dynamics.  Even 
an  acoustic  wave  at  very  low  amplitude  could  lead  to  noticeable  sinuous  wave-like  structures  on 
the  jet  surface.  The  two-dimensional  characteristics  of  those  wavy  structures  are  confirmed. 
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Task  4 

Flow  and  Flame  Dynamics  of  Shear  Co-axial  Injectors  with 
Liquid  Oxygen  (LOX)  and  Methane 

The  numerical  analysis  has  also  been  applied  to  examine  the  key  phenomena  and 
mechanisms  associated  with  the  coaxial  injection  and  combustion  of  LOX  and  methane.  Emphasis 
is  placed  on  the  flame  stabilization  and  flow/flame  interactions  in  the  near  field.  Co-flowing 
methane  (outer)  and  oxygen  (inner)  streams,  separated  by  the  central  LOX  post,  are  injected 
through  the  injector.  The  inner  diameters  of  the  LOX  post  and  the  methane  annulus  are  3.4  and  5.2 
mm,  respectively.  Fully  developed  turbulent  pipe  flows  are  assumed  at  the  injector  exit.  The 
computational  domain  downstream  of  the  injector  measures  a  length  of  40 DL0X  and  a  radius  of 

12 Dwx ,  where  Dmx  is  the  diameter  of  the  LOX  tube.  Such  dimensions  are  sufficient  to  minimize 
the  effect  of  the  far-field  boundary  conditions  on  the  near-field  flow  evolution.  A  quasi- 
axisymmetric  simulation  was  conducted.  The  entire  grid  system  consists  of  360x240  cells  along 
the  axial  and  radial  directions,  respectively.  The  mean  grid  size  of  10  pm  falls  in  the  subinertial- 
range  of  the  turbulent  kinetic  energy  spectrum  estimated  based  on  the  inlet  Reynolds  number  of  the 
oxygen  stream.  To  resolve  the  wake  region  and  vortical  structures  which  shed  downstream  from 
the  rim  of  the  LOX  post,  a  total  of  40  grid  points  are  used  to  cover  the  LOX  post  in  the  radial 
direction.  Two  different  simulation  conditions  are  considered  herein.  In  Case  1,  liquid  oxygen  and 
methane  are  injected  at  temperatures  of  122  and  300  K,  respectively.  The  bulk  velocities  of  the  two 
streams  are  13  and  75  m/s,  respectively.  The  mixture  ratio  of  3  is  typical  of  operational  engines.  In 
Case  2,  the  LOX  inlet  temperature  and  velocity  remain  fixed,  but  the  temperature  of  the  methane 
stream  increases  to  450  K.  Consequently,  the  methane  stream  velocity  becomes  132  m/s  to 
maintain  the  same  mixture  ratio  as  that  in  Case  1.  The  momentum  flux  ratio  of  the  two  streams, 
increases  from  2.5  of  Case  1  to  4.3  of  Case  2.  Such  a  variation  allows  a  careful  investigation  of  the 
effects  of  momentum-flux  ratio,  a  key  parameter  in  the  design  of  shear  coaxial  injectors,  on  the 
flame  stabilization  characteristics. 

Results  from  this  task  have  led  to  the  following  technical  publications. 

“Near-Field  Flow  and  Flame  Dynamics  of  LOX/Methane  Shear  coaxial  Injector  under 

Supercritical  Conditions,”  by  N.  Zong  and  V.  Yang,  to  appear  in  Proceedings  of  the 

Combustion  Institute,  2006. 

“Supercritical  LOX/Methane  Flame  Stabilization  and  Dynamics  in  a  Shear  Coaxial 

Injector,”  by  N.  Zong  and  V.  Yang,  AIAA  Paper  2006-0760,  presented  at  the  44th  AIAA 

Aerospace  Sciences  Meeting  and  Exhibit,  January  2006. 

“High-  Pressure  LOX/Methane  Mixing  and  Combustion  Processes,”  by  N.  Zong  and  V. 

Yang,  AIAA  Paper  2005-0152,  presented  at  the  43rd  AIAA  Aerospace  Sciences  Meeting 

and  Exhibit,  January  2005. 

4.1  Introduction 

Extensive  experimental  and  theoretical  studies  have  been  conducted  to  improve  the 
understanding  of  cryogenic-propellant  mixing  and  combustion  under  both  trans-  and  super-critical 
conditions.  Most  of  the  previous  work  was  focused  on  systems  involving  liquid  oxygen  (LOX)  and 
gaseous  hydrogen  [1].  For  future  development  of  high-performance  reusable  liquid  rocket  engines, 
combustion  of  LOX  and  methane  has  recently  attracted  considerable  interest.  The  key  issues  of 
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concern  include  injector  dynamics,  combustion  efficiency  and  stability,  and  chamber  cooling,  to 
cite  a  few. 

As  an  attempt  to  address  these  questions,  Zurback  et  al.  [2]  conducted  a  preliminary  flow 
visualization  of  shear  coaxial  injection  and  combustion  of  LOX  and  methane  at  near-critical 
pressures.  Shadowgraph  images  were  obtained  to  reveal  that  the  flame  was  attached  to  the  LOX 
post  and  spreaded  downstream  along  the  oxygen  jet  boundary,  a  phenomenon  quite  similar  to  that 
observed  for  the  LOX/hydrogen  system.  Singla  et  al.  [3]  later  performed  a  more  detailed 
experiment  of  high-pressure  oxygen  and  methane  combustion  associated  with  a  shear  coaxial 
injector.  The  temperature  of  the  oxygen  stream  remained  at  85  K,  whereas  the  methane  stream  took 
a  value  between  120  and  288  K  to  simulate  trans-  and  super-critical  injection  conditions.  The 
chamber  pressure  varied  from  4.5  to  6.0  MPa.  Emission  images  of  excited  OH  and  CH  radicals 
were  recorded  and  time-averaged  to  determine  the  mean  flame  structure.  Results  indicated  that  the 
flame  was  stabilized  in  the  vicinity  of  the  LOX  post  tip  under  all  flow  conditions.  Since 
vaporization  was  the  slowest  process  at  a  subcritical  temperature,  part  of  the  unbumed  oxygen 
droplets  penetrated  into  the  inner  flame.  After  vaporization  and  mixing  with  gaseous  methane  at 
the  outer  boundary  of  the  methane  stream,  a  second  flame  with  a  greater  expansion  angle  was 
formed.  Therefore,  when  both  LOX  and  methane  were  injected  at  a  transcritical  condition,  the 
flame  featured  two  different  regions  of  light  emission,  one  surrounding  the  liquid  oxygen  jet  and 
the  other  located  close  to  the  outer  boundary  of  the  annular  methane  stream.  The  situation  changed 
if  oxygen  was  injected  at  a  subcritical  temperature  while  methane  was  gaseous.  The  resultant 
enhancement  of  turbulent  mixing  led  to  only  one  flame  surrounding  the  oxygen  jet. 

In  parallel  to  the  experimental  studies,  a  comprehensive  numerical  model  for  high-pressure 
fluid  mixing  and  combustion  was  recently  developed  [4],  The  analysis  is  employed  in  the  present 
work  to  explore  various  fundamental  physiochemical  mechanisms  associated  with  shear  co-axial 
mixing  and  combustion  of  LOX  and  methane  under  representative  engine  operating  conditions. 
The  effects  of  the  momentum  flux  ratio  of  the  two  propellant  streams  on  the  near-injector  flow  and 
flame  dynamics  are  examined  in  detail. 

4.2  Theoretical  formulation  and  numerical  treatment 

The  basis  of  the  present  work  is  the  general  theoretical/numerical  framework  described  in 
Refs.  [5-678].  The  formulation  accommodates  the  full  conservation  laws  and  real-fluid 
thermodynamics  and  transport  over  the  entire  temperature  and  pressure  regime  of  concern. 
Turbulence  closure  is  achieved  by  means  of  a  large-eddy-simulation  (LES)  technique,  in  which 
large-scale  motions  are  calculated  explicitly  and  the  effects  of  unresolved  small-scale  turbulence 
are  modeled  either  analytically  or  empirically.  The  Favre-filtered  mass,  momentum,  energy,  and 
species  conservation  equations  are  derived  by  filtering  small-scale  dynamics  from  resolved  scales 
over  a  well-defined  set  of  spatial  and  temporal  intervals.  The  effects  of  subgrid-scale  (sgs)  motions 
are  treated  using  the  model  proposed  by  Erlebacher  et  al.  [9].  It  employs  a  Favre-averaged 
generalization  of  the  Smagorinsky  eddy  viscosity  model  coupled  with  a  gradient-diffusion 
assumption  to  simulate  sgs  energy  and  species  transport  processes.  The  Smagorinsky  coefficients 
CR  (=0.01)  and  Cl  (=0.007)  are  determined  empirically.  Thermodynamic  properties,  such  as 
enthalpy,  Gibbs  energy,  and  constant-pressure  specific  heat,  are  obtained  directly  from  fundamental 
thermodynamics  theories  and  a  modified  Soave-Redlich-Kwong  (SRK)  equation  of  state  [10]. 
Transport  properties,  such  as  viscosity  and  thermal  conductivity,  are  evaluated  using  an  extended 
corresponding-state  theory  [11,  12]  along  with  the  32-term  Benedict-Webb-Robin  equation  of  state 
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[13].  Mass  diffusivity  is  obtained  by  means  of  the  Takahashi  method,  calibrated  for  high  pressure 
conditions  [14].  The  implementation  and  validation  of  the  property  evaluation  schemes  were 
detailed  by  Yang  [15]  and  Meng  et  al.  [8].  Owing  to  the  lack  of  detailed  chemical  kinetics  schemes 
for  oxygen  and  methane  at  high  pressures,  a  one-step  global  reaction  model  involving  four  species 
(CH4,  02,  C02,  H20)  is  employed  for  the  combustion  of  oxygen  and  methane  [16]. 

4.2.1  Turbulence/flame-structure  interaction 

The  modeling  of  turbulence/chemistry  interaction  in  a  physically  meaningful  manner 
represents  a  critical  and  challenging  issue  in  the  present  study  of  supercritical  combustion.  Cuenot 
and  Poinsot  [17]  suggested,  based  on  the  results  from  the  direct  numerical  simulations  (DNS)  of 
flame/vortex  interactions,  that  such  interaction  in  a  non-premixed  turbulent  flame  be  characterized 
by  two  non-dimensional  parameters:  turbulent  Reynolds  number  Re  ,  and  turbulent  Damkohler 

number  n  ,  defined  respectively  as  follows. 

Re,  =  ——  0  (/,  /  ld  )4/3  (1) 


A,  =-  =  --□  2^Af 

Tc  Tc  (2) 

where  rf,  tc,  and  rk  are  the  characteristic  times  for  turbulent  integral-scale  motions,  chemical 
reactions,  and  molecular  diffusion,  respectively,  and  v  the  kinematic  viscosity.  The  turbulent 
integral  and  Kolmogorov  length  scales  are  denotes  by  /,  and  ld,  respectively,  u  is  the 

characteristic  velocity  of  eddies  of  size  /, .  The  local  flame  Damkohler  number  Dd  is  defined  as 
the  ratio  of  the  flow  residence  time  to  the  chemical  time  ( Dd  -  tflTc ).  Depending  on  the  values  of 

the  Reynolds  and  Damkohler  numbers,  the  characteristics  of  non-premixed  combustion  can  be 
classified  into  four  different  regimes  as  illustrated  in  Fig.  1  [17].  Thin  flames,  also  identified  as 
flamelets,  take  place  when  the  chemical  time  tc  is  small  and  Damkohler  number  Da  is  large.  As 

the  chemical  time  tc  increases  and  the  flame  thickness  becomes  the  same  order  of  magnitude  as  the 

length  scale  of  the  Kolmogorov  eddies,  the  emerging  of  unsteady  effects  must  be  taken  into 
account.  Flame  extinction  occurs  for  very  slow  chemical  reactions,  corresponding  to  the  low 
Damkohler  number  regime  in  Fig.  1.  For  a  turbulent  diffusion  flame,  the  local  flame  thickness  and 
reaction  rate  depend  strongly  on  such  flow  conditions  as  the  strain  rate,  and  are  affected  by  various 
unsteady  effects  encountered.  The  flame  structures  may  vary  at  different  spatial  locations  in  the 
flowfield. 
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Fig.  1  Regimes  for  turbulent  non-premixed  combustion  as  a  function  of  the  Damkohler  number 
Da—Tt  I  Tc  and  the  turbulent  Reynolds  number  Re( . 


Several  closure  schemes,  such  as  the  linear-eddy  [18],  conditional-momentum-closure  [19], 
laminar-flamelet  [20],  and  probability-density-function  (PDF)  [21]  methods,  have  been  proposed. 
An  accurate  and  efficient  treatment  of  turbulence/chemistry  interactions  within  the  context  of  LES, 
however,  still  remains  to  be  established,  even  for  ideal-gas  mixtures.  The  situation  becomes  more 
challenging  for  supercritical  fluids  due  to  complications  arising  from  steep  property  variations  and 
the  lack  of  basic  flame  properties  at  high  pressures.  Even  a  well-calibrated  sgs  model  for  pure 
fluid-dynamic  processes  is  not  currently  available.  In  light  of  this  limitation,  a  simplified  approach 
is  adopted  in  the  present  work.  It  is  well  known  that  the  Reynolds  number  increases  almost  linearly 
with  pressure.  An  increase  in  pressure  from  1  to  100  atm  gives  rise  to  an  increase  of  the  Reynolds 
number  by  two  orders  of  magnitude.  The  corresponding  Kolmogorov  microscale  increases  by  1.5 
orders  of  magnitude.  As  a  consequence,  under  supercritical  conditions,  the  Reynolds  number  may 
reach  such  a  level  that  the  turbulent  eddies  may  penetrate  into  the  flame  zone  and  greatly  enhance 
the  mixing  process.  This  leaves  chemical  reactions  a  rate-controlling  process.  The  resolved-scale 
chemical  source  can  thus  be  evaluated  through  a  direct  closure  approach  without  considering  the 
contributions  from  sgs  fluctuations.  The  filtered  reaction  rate  co  is  modeled  as. 


Qjic(p,T,YlJ2,-JN)  =  d\(p,T,Yi,Y2, -,YN) 


(3) 


In  spite  of  the  neglect  of  sgs  contributions,  the  approach  allows  for  fluctuations  in  species 
production  rates  as  functions  of  instantaneous  flow  properties  and  accounts  for  finite-rate 
chemistry.  The  model  uncertainty  and  error  can  be  effectively  reduced  if  the  spatial  and  temporal 
resolutions  in  numerical  simulations  approach  the  length  and  time  scales  of  sgs  fluctuations. 


4.2.2  Computational  domain  and  boundary  conditions 

Figure  2  shows  the  physical  model  under  consideration.  A  methane  (outer)  and  an  oxygen 
(inner)  stream,  separated  by  a  0.38  mm  thick  LOX  post,  are  delivered  to  a  co-axial  injector.  The 
inner  diameter  of  the  LOX  post  is  3.42  mm,  and  the  outer  diameter  of  the  methane  annulus  is  5.18 
mm.  The  injector  geometry  is  chosen  to  match  that  employed  in  an  experimental  study  of  high- 
pressure  LOX/methane  combustion  [22].  The  computation  domain  includes  the  injector  and  a 
downstream  region  that  measures  6 Dwx  and  40-D^  in  the  radial  and  axial  directions, 
respectively.  Because  of  the  enormous  computational  effort  required  for  calculating  the  flowfield  in 
the  entire  three-dimensional  regime,  only  a  cylindrical  sector  with  periodic  boundary  conditions 
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specified  in  the  azimuthal  direction  is  treated  herein.  The  analysis,  in  spite  of  the  lack  of  vortex¬ 
stretching  mechanism,  has  been  shown  to  be  able  to  capture  the  salient  features  of  supercritical 
fluid  injection  and  mixing  dynamics  [4]. 
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*  rL 
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i- 


DCH4=5.18  mm 


40D 
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Fig.  2  Schematic  diagram  of  a  shear  coaxial  injector. 


The  velocity  profile  for  a  fully  developed  turbulent  pipe  flow  is  assumed  at  the  injector  inlet 
for  both  streams.  Turbulence  is  provided  by  superimposing  broad-band  white  noise  onto  the  mean 
velocity  profile.  The  disturbances  are  generated  by  a  Gaussian  random-number  generator  with  an 
intensity  of  1%  of  the  mean  quantity.  The  inflow  perturbations  do  not  contain  any  dominant 
frequency  and  only  serve  to  trigger  the  flow  instabilities  inherent  in  the  shear  layers.  At  the 
downstream  boundary,  extrapolation  of  primitive  variables  from  the  interior  may  cause  undesired 
reflection  of  waves  propagating  into  the  computational  domain.  Thus,  the  non-reflecting  boundary 
conditions  based  on  the  method  of  characteristics  are  applied,  along  with  the  specification  of  a 
reference  pressure  [4],  At  the  radial  boundaries,  the  pressure  and  temperature  are  specified  as  the 
ambient  values.  The  axial,  radial,  and  azimuthal  velocities  are  extrapolated  from  the  interior. 
Finally,  the  no-slip  adiabatic  condition  is  enforced  along  the  solid  walls. 

4.2.3  Numerical  framework 

The  theoretical  formulation  outlined  above  requires  a  robust  computational  scheme,  due  to 
the  numerical  stiffness  caused  by  rapid  flow  property  variations  and  wide  disparities  of  the 
characteristic  time  and  length  scales  involved.  To  this  end,  a  unified  treatment  of  general  fluid 
thermodynamics,  based  on  the  concepts  of  partial-mass  and  partial-density  properties,  and  valid  for 
the  entire  pressure  and  temperature  regimes  of  concern,  is  established  and  incorporated  into  a 
preconditioning  scheme  [6].  Since  the  numerical  relations,  including  the  Jacobian  matrices  and 
eigenvalues,  are  derived  directly  from  fundamental  thermodynamics  theories,  the  algorithm  is  self- 
consistent  and  robust. 

The  numerical  framework  employs  a  density-based,  finite  volume  methodology  along  with 
a  dual  time-step  integration  technique  [23],  Temporal  discretization  is  obtained  using  a  second- 
order  backward  difference,  and  the  inner-loop  pseudo-time  term  is  integrated  with  a  four-step 
Runge-Kutta  scheme.  Spatial  discretization  is  achieved  with  a  fourth-order,  central-difference 
scheme  in  generalized  coordinates.  In  the  present  study  of  LOX/methane  combustion,  exceedingly 
large  property  gradients  exist  in  the  near  field  of  the  injector  where  the  fluid  density  varies  by  a 
factor  of  1000  over  a  fraction  of  the  LOX  post  thickness  (0.38  mm).  To  effectively  suppress 


74 


numerical  oscillations  in  this  region  and  ensure  computational  stability,  the  fourth-order  scalar 
dissipation  with  a  total-variation-diminishing  switch  developed  by  Swanson  and  Turkel  [24]  is 
implemented.  The  dissipation  coefficient  is  selected  to  be  £4  =0.001  to  minimize  the 
contamination  from  numerical  dissipation. 

4.3  Results  and  discussions 

The  theoretical  model  and  numerical  scheme  established  in  the  preceding  sections  are 
implemented  to  study  the  coaxial  injection  and  combustion  of  LOX  and  methane  under  supercritical 
conditions.  Table  1  summarizes  the  flow  conditions  considered  herein.  In  Case  1,  liquid  oxygen 
and  methane  are  injected  at  temperatures  of  122  and  300  K,  respectively.  The  bulk  velocities  of  the 
two  streams  are  13  and  75  m/s,  respectively.  The  mixture  ratio  of  3  is  typical  of  operational 
engines.  In  Case  2,  the  LOX  inlet  temperature  and  velocity  remain  fixed,  but  the  temperature  of  the 
methane  stream  increases  to  450  K.  Consequently,  the  methane  stream  velocity  becomes  132  m/s 
to  maintain  the  same  mixture  ratio  as  that  in  Case  1.  The  momentum  flux  ratio  of  the  two  streams, 
thus,  increases  from  2.5  of  Case  1  to  4.3  of  Case  2.  Such  a  variation  allows  a  careful  investigation 
of  the  effects  of  momentum-flux  ratio,  a  key  parameter  in  the  design  of  shear  coaxial  injectors  [25], 
on  the  flame  stabilization  characteristics.  It  is  worth  mentioning  that  for  both  cases,  the  methane 
stream  remains  at  a  supercritical  state,  while  the  LOX  stream  enters  the  injector  at  a  transcritical 
temperature.  The  volume  downstream  of  the  injector  is  preconditioned  with  H20  and  C02  with  a 
stoichiometric  composition  for  oxygen/methane  combustion.  The  chamber  pressure  of  10  MPa  is 
well  above  the  critical  pressures  of  oxygen  (5.04  MPa)  and  methane  (4.60  MPa). 


Table  1  Simulation  Conditions 


Case  1 

Case  2 

p  (atm) 

100 

100 

Tlox  (K) 

122 

122 

Tch,  (K) 

300 

450 

Pwx  (kg/m3) 

1006 

1006 

Pcna  (kg/m3) 

75 

42 

^LOX  (m/s) 

13.06 

13.06 

uCHs  (m/s) 

75 

132 

(Pu  )ch4  f(Pu  )lox 

2.5 

4.3 

™LOX  CHa 

3.0 

3.0 

Rclox 

3.9xl05 

3.9x1 0s 

ReCHA 

1.2x10s 

9.5  xlO4 

The  computational  grid  for  the  regime  downstream  of  the  injector  consists  360x240  points 
along  the  axial  and  radial  directions,  respectively.  The  grids  are  clustered  in  the  shear  layers  and 
near  the  injector  to  resolve  rapid  property  variations  in  those  regions.  The  smallest  grid  size  in  the 
radial  direction  is  6  pm ,  which  well  falls  in  the  inertial  sub-range  of  the  turbulent  kinetic  energy 
spectrum  estimated  using  the  Kolmogorov-Obukhow  theory.  The  computational  domain  is  divided 
into  31  blocks,  with  each  calculated  on  a  single  processor  of  a  distributed  computing  facility.  The 
physical  time  step  is  100  nano-seconds  and  the  maximum  CFL  number  for  the  inner-loop  pseudo- 


75 


time  integration  is  0.7.  A  grid  independence  study  is  performed  on  a  fine  grid  of  540x360  cells 
for  Case  1.  No  discernible  changes  of  the  near-field  flow  and  flame  dynamics  are  observed  in 
terms  of  the  mean  properties  and  vortex  shedding  frequency  phenomenon.  The  grid  system 
employed  in  the  present  work  appears  to  be  credible. 

Figure  3  shows  the  instantaneous  fields  of  temperature,  density,  vorticity,  and  mass  fractions 
of  methane,  oxygen,  and  water  in  the  vicinity  of  the  LOX  post.  A  diffusion-dominated  flame 
emanates  immediately  from  the  LOX  post  and  propagates  downstream  along  the  surface  of  the 
LOX  stream.  A  wake  region,  which  consists  of  hot  combustion  products,  effectively  separates  the 
methane  and  LOX  streams.  Similar  to  the  non-reacting  flow  cases  [26],  the  near-field  flow 
dynamics  are  characterized  by  the  evolution  of  the  three  mixing  layers  originating  from  the  inner 
and  outer  edges  of  the  methane  annulus  and  the  inner  rim  of  the  LOX  post.  The  development  of  the 
inner  mixing  layer  of  the  methane  stream,  however,  is  slightly  inhibited  by  the  expansion  of  the 
combustion  products  in  the  flame  zone.  Because  of  the  liquid-like  behavior  of  the  oxygen  jet,  a 
steep  density  gradient  exists  between  the  flame  and  the  oxygen  stream.  As  a  result,  the  large-scale 
vortices  emerging  from  the  outer  rim  of  the  LOX  post  evolve  in  a  manner  analogous  to  that 
produced  by  a  backward-facing  step  and  mainly  reside  on  the  lighter  fluid  side.  The  evolutions  of 
those  vortices  significantly  enhance  the  mixing  of  methane  and  hot  products.  The  denser  oxygen 
stream  is  less  influenced.  The  vortices  interact  and  coalescence  with  their  neighboring  vortices 
while  convecting  downstream.  The  flame  residing  between  the  oxygen  and  methane  streams 
appears  to  be  quite  resistant  to  flow  straining.  Even  the  strong  strain  rate  generated  during  the 
vortex -paring  process  does  not  cause  the  flame  front  to  break  up. 

The  higher  momentum-flux  ratio  in  Case  2  exerts  a  substantial  influence  on  the  flow 
evolution.  Both  the  inner  and  outer  mixing  layers  of  the  methane  stream  become  more  dynamic, 
and  the  associated  instability  waves  roll  up  into  vortices  at  an  upstream  location.  The  vortices 
formed  downstream  of  the  LOX  post  are  stronger,  thereby  enhancing  the  mixing  between  methane 
and  hot  combustion  products  and  shortening  the  potential  core  of  the  LOX  stream.  The  interactions 
among  vortices  give  rise  to  a  much  more  complicated  flowfield  with  the  emergence  of  small 
structures.  The  potential  core  of  the  methane  stream  is  also  significantly  reduced. 

To  further  explore  the  flow  structure,  the  instantaneous  flow  properties  are  time  averaged 
over  two  flow-through  times  after  the  calculated  flowfield  has  reached  its  stationary  state.  Figure  4 
shows  the  time-averaged  fields  of  temperature,  density,  specific  heat,  and  compressibility  factor  of 
Case  1.  The  constant-pressure  specific  heat  has  been  normalized  with  respect  to  the  value  at  the 
initial  state  of  the  LOX  jet.  The  potential  core  of  the  methane  stream  and  the  flame  zone  are  clearly 
observed  in  the  temperature  field.  The  fluid  density  varies  from  10  to  1000  kg/m3  over  a  very 
small  spatial  domain  close  to  the  LOX  post.  Such  a  steep  density  gradient  makes  this  region 
behave  like  a  contact  discontinuity  [5,27].  The  specific  heat  exhibits  an  anomalous  variation  across 
the  LOX  jet  boundary.  It  increases  rapidly  and  reaches  a  peak  as  the  fluid  temperature  transits  from 
a  subcritical  to  a  supercritical  value  on  the  isobaric  curve  [7].  This  effect,  combined  with  the  low 
thermal  conductivity  (not  shown),  facilitates  the  formation  of  a  strong  density-gradient  region 
between  the  LOX  jet  and  hot  combustion  products  [7].  The  compressibility  factor,  which  measures 
the  departure  from  the  ideal-gas  behavior,  has  values  of  0.30  and  0.87  in  the  core  regions  of  the 
LOX  and  methane  streams,  respectively.  It  increases  expeditiously  and  approaches  unity  in  the  hot 
products,  where  the  real-fluid  effects  essentially  diminish  due  to  the  local  high  temperature. 


Case  1 


Case  2 
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Fig.  3  Snapshots  of  distributions  of  temperature,  density,  vorticity,  and  mass  fractions  of  oxygen, 
methane,  and  water  immediate  downstream  of  the  injector  LOX  post  for  Cases  1  and  2  in  Table  1 


Fig.  4  Short-time  averaged  temperature,  density,  constant-pressure  specific  heat,  and 
compressibility-factor  fields  of  Case  1  in  Table  1. 


Flame  stabilization  represents  a  critical  issue  in  the  combustor  design.  For  a  shear  coaxial 
injector  with  cryogenic  propellants,  the  flame  is  stabilized  by  the  recirculation  flow  downstream  of 
the  LOX  post,  a  phenomenon  first  quantified  in  Ref.  [5].  Several  experimental  and  numerical 
studies  [5,  27-29]  of  shear  coaxial  injection  and  combustion  of  LOX/hydrogen  have  confirmed  that 
such  a  strong  recirculating  flow  acts  as  a  hot-product  pool  providing  the  energy  to  ignite  incoming 
propellants.  Figure  5  shows  the  temporal  evolution  of  the  temperature  and  velocity-vector  fields  in 
the  vicinity  of  the  LOX  post.  The  solid  lines  in  the  velocity-vector  fields  correspond  to  the 
isotherm  of  3000  K,  which  represents  the  flame  boundary.  The  large-scale  vortices  emerging  from 
the  outer  rim  of  the  LOX  post  facilitate  the  mixing  between  the  incoming  methane  stream  and  hot 
products.  Driven  by  those  vortices,  a  relative  weak  recirculation  flow  forms  close  to  the  inner  rim 
the  LOX  post  and  convects  the  oxygen-rich  products  toward  the  methane  stream.  Unlike  the  case 
with  LOX/hydrogen  combustion,  in  which  the  flame  is  anchored  very  close  to  the  LOX  jet 
boundary  due  to  the  high  diffusivity  of  hydrogen  and  the  strong  inertia  of  the  LOX  jet  [5,27],  the 
LOX/methane  flame  is  anchored  between  two  counter-rotating  wake  recirculation  zones,  as 
illustrated  schematically  in  Fig.  6. 


Fig.  5  Time  evolution  of  temperature  and  velocity-vector  fields  in  the  vicinity  of  LOX  post  of  Case 
1  in  Table  1;  solid  lines  correspond  to  the  isotherm  of  3000  K. 


methane 


Fig.  6  Schematic  diagram  of  flame  anchoring  mechanism. 
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Figure  7  shows  the  time  histories  of  the  fluctuating  axial  velocities  at  three  different  axial 
locations  along  the  inner  and  outer  shear  layers  of  the  methane  stream.  The  locations  of  Probes  1, 
2,  and  3  correspond  to  the  positions  of  the  initial  instability  wave,  vortex  roll-up,  and  vortex  pairing 
in  the  outer  shear  layer.  Small  velocity  fluctuations  are  observed  in  the  early  stage  of  the  shear 
layer  development.  As  the  vortices  roll-up  downstream,  periodic  flow  oscillations  are  clearly 
observed  at  Probe  2.  The  amplitudes  of  the  velocity  fluctuations  further  increase  with  the  growth  of 
the  vortices  (Probe  3).  The  effects  of  chemical  reactions  on  the  vortex  evolution  can  be  identified 
from  the  temporal  behaviors  at  Probes  4-6  along  the  inner  shear  layer,  which  are  close  to  the  flame 
zone.  In  addition  to  those  well-organized  flow  oscillations,  there  exist  high-frequency  fluctuations. 
The  phenomenon  can  be  attributed  to  the  volume  dilatation  produced  by  the  oscillatory  flame  [30]. 
As  a  consequence  of  the  strong  interactions  between  vortices,  as  well  as  their  coupling  with  the 
flame,  the  velocity  fluctuation  at  Probe  6  starts  to  loss  its  periodic  identity.  A  closely  similar 
behavior  is  observed  in  Case  2  (not  shown). 
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Fig.  7  Time  histories  of  axial  velocity  oscillations  along  the  inner  and  outer  mixing  layers  of  the 
methane  stream  of  Case  1  in  Table  1. 

Figure  8  shows  the  power  spectral  densities  of  the  axial-velocity  fluctuations  at  Probe  2  and  5 
for  Cases  1  and  2.  A  dominant  frequency  of  13.8  KHz,  corresponding  to  the  vortex  shedding 
frequency,  is  obtained  in  both  the  outer  and  inner  shear  layers  for  Case  1.  Higher  harmonics  exist 
in  the  outer  shear  layer,  but  their  magnitudes  become  much  smaller  in  the  inner  shear  layer.  For 
Case  2,  the  higher  momentum  flux  of  the  methane  stream  causes  an  increase  of  the  vortex  shedding 
frequency  to  17.2  kHz  .  Rehab  et  al.  [  31]  studied  the  flow  characteristics  of  a  coaxial  water  jet 
without  a  splitter  between  the  two  streams.  They  suggested  that  the  outer  shear-layer  dynamics  are 
dominant  over  the  inner  one,  and  the  most  amplified  frequency  can  be  predicted  by  the  correlation 
proposed  by  Schadow  and  Gutmark  [32],  The  vortex  shedding  frequencies  obtained  in  the  present 
work,  however,  do  not  confirm  their  correlation.  The  existence  of  a  second  length  scale  (i.e.,  the 
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Case  1 


Case  2 


Fig.  8  Frequency  spectra  of  axial  velocity  oscillations  in  the  inner  and  outer  mixing  layers  of 
the  methane  stream  of  Cases  1  and  2  in  Table  1. 

LOX  post  thickness)  and  the  severe  density  stratification  between  the  LOX  jet  and  hot  products 
significantly  modify  the  near-field  flow  dynamics.  The  vortices  shed  from  the  LOX  post  tip  in  a 
manner  similar  to  the  vortex  generated  behind  a  backward-facing  step,  where  the  eddy-formation 
frequency  is  on  the  order  of  0.1  in  terms  of  the  Strouhal  number  defined  based  on  the  bulk  velocity 
of  the  outer  stream  and  the  step  thickness  [33]. 

4.4  Conclusions 

A  comprehensive  theoretical/numerical  framework  has  been  established  to  study  the  mixing 
and  combustion  of  co-flowing  liquid  oxygen  (LOX)  and  gaseous  methane  through  a  shear  co-axial 
injector.  The  physical  model  and  flow  conditions  are  representative  of  cryogenic-propellant  rocket 
engine  injectors  operating  at  supercritical  pressures.  The  formulation  is  based  on  the  full 
conservation  laws,  and  takes  into  account  real-fluid  thermodynamics  and  transport  phenomena. 
Turbulence  closure  is  achieved  using  a  large-eddy-simulation  technique.  The  work  allows  for  a 
thorough  investigation  into  the  injector  flow  development  and  flame  dynamics.  As  a  consequence 
of  the  large  density  stratification,  the  flame  is  anchored  in  the  wake  recirculating  region  behind  the 
LOX  post  and  propagates  along  the  boundary  of  the  LOX  jet.  The  near-field  flow  evolution  is 
dictated  by  the  large-scale  vortices  in  the  inner  shear  layer  of  the  methane  stream.  The  dominant 
frequencies  of  vortex  shedding  match  those  for  a  rear-facing  step  flow,  a  phenomenon  that  can  be 
attributed  to  the  strong  inertia  of  the  LOX  stream  and  lighter  density  of  gaseous  methane.  The 
effects  of  the  momentum-flux  ratio  of  the  two  streams  on  the  injector  flow  and  flame  characteristics 
are  also  examined. 
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